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Abstract 



The muon decay parameter delta characterizes momentum dependence of the parity- 
violating muon decay asymmetry. A new measurement of delta has been performed 
using the first physics data recorded by the TWIST experiment at TRIUMF. The 
obtained value, 5 = 0.74964 ±0.00066 (stat.) ±0.00112 (syst.), is consistent with the 
Standard Model expectation 5 = 3/4. This is the first determination of 5 performed 
using a blind analysis technique. Combined with other data, the measurement 
sets new model-independent limits on effective right-handed couplings of the muon. 
Improved limits on the product of another muon decay parameter, £, and the muon 
polarization in pion decay, P M , are obtained in the form: 0.9960 < < £ < 1.0040, 
at 90% confidence level. Implications for left-right symmetric models are discussed. 
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Chapter 1 

Introduction 



The Standard Model of particle physics [H [3 [3] has been very successful in de- 
scribing available experimental data. The only observed deviation from the original 
theory is the recent observation of neutrino oscillations, but that can be easily ac- 
commodated in the model's framework and does not lead to a conceptual change. 
Despite the many successes, the theory is generally believed to be incomplete, see 
e.g. Section VII in pQ. Numerous extensions of the Standard Model have been 
proposed, and experimental searches for New Physics are ongoing. The experimen- 
tal efforts explore two complementary approaches. One type of experiment aims at 
direct observation of new particles. These searches require the energy of the col- 
lision to be high enough to produce the supposed heavy particle being looked for, 
and their reach is limited by the capabilities of the accelerator. The other direction 
of such research exploits contributions of hypothetical particles to known processes 
through virtual (loop) effects. These experiments can be done at low energies. The 
mass-scale reach of this kind of search is limited by the precision of the measurement 
and by the theoretical precision of the calculation of the "known" processes. 

Muon decay \x — > evv, studied by TWIST (TRIUMF Weak Interaction Symme- 
try Test experiment), is one of a few processes in particle physics that can be un- 
ambiguously calculated with high accuracy in the framework of a theoretical model. 
The purely leptonic nature of the decay eliminates many uncertainties due to the 
internal structure of the particles. The strong interaction, which at present can not 
be accurately evaluated from first principles, enters only through higher order radia- 
tive corrections. The fractional hadronic contribution to the energy spectrum can 
be estimated as 0.07 (a/vr) 2 « 0.4 x 10" 6 [5], so any related uncertainty is negligible 
for the current state of the field. On the other hand muons are easy to produce 
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in large quantities at an accelerator. That means high experimental statistics is 
affordable, and precision experiments can be done to test theoretical predictions. 
Thus the decay of a muon is an ideal low energy process with which to investigate 
the Lorentz structure of the weak interaction. 

1.1 The 4-fermion interaction formalism 

An approach useful in searches for new physics is to start with a very general de- 
scription of the process, then try to limit the possibilities. The most general, local, 
derivative- free, Lorentz-invariant and lepton-number conserving four fermion inter- 
action was introduced by Michel [6]. Using the notation of [7J, which represents 
particles by fields of definite chirality [HJ [9j [10] , the interaction matrix element can 
be written in a "helicity projection form" as 

M = ^Jf £ ^(e e |n|(^ e ) n )((^) m |r 7 |^). (1.1) 
1 = S,V,T 
e, fj, = R, L 

Here Gf is the Fermi coupling constant, while 7 labels scalar, vector or tensor type 
of interaction: 

V s = 1, T v = 7 a , T T = ^ r a^. 

a/2 

In the last equation 7° are the Dirac gamma matrices, and a a ^ = §(7 Q 7^ — 7^7°) • 
The indices e and /i indicate the chirality (handedness) of the spinors of the charged 
leptons: 

^R,L = X (1 ± 75) i> 

The chiralities n and m of the u e and spinors, respectively, are uniquely deter- 
mined for given 7, e and fi. The tensor term in (jl.ip requires special attention. Due 
to the identity 75 a a p = | e a p\ p a Xp , the coupling constants gJ> R = g\ L = 0. So 
the general interaction (jl.ip is defined by 10 complex parameters. Since a common 
phase does not matter, the interaction is fully described by 19 real independent 
coupling constants. The usual convention is to absorb the overall strength of the 
interaction into Gf, and normalize gl^ |llj as: 

ll^l 2 + Ml? + Mr\ 2 + Ml? 

+ Kr? + \9 V RL? + \9 V LR? + \9 V LL? (1-2) 
+ 3|^ L | 2 + Z\g T LR \ 2 = 1. 
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Using (jl.ip as a starting point, it is straightforward to calculate the differential 
rate, in energy and angle, of positrons emitted in muon decay \12 \ 113 } fT4"]: 

dx f c i {e) =^W^G%^4{F ls (x) + P,cos(e)F AS ( X )} (1.3) 

with 

Fjs(x) = x (1 - x) + -p (Ax 2 - 3x - xg) + 7?x (1 - x) + Fi| C (x), (1.4) 

Fas(x) = ^i\jx 2 - xl 

(1.5) 

Here is the muon mass, W efl = [rr? + m 2 )/2m^ is the maximum energy of 
the emitted positron, x = E e /W e ^ is the reduced positron energy, is the angle 
between the positron momentum and an arbitrary direction z, xq = m e /W e ^ is the 
dimensionless electron mass, — 1 < < 1 is the muon polarization with respect to 
z, and Fjg (x) and F\s{x) are radiative corrections. The muon decay parameters 
p, 7], £, S, are real numbers expressed through bilinear combinations of the coupling 
constants g]^, and the indices IS and AS label the isotropic and anisotropic terms. 

In the Standard Model muon decay is mediated by a W vector boson. It is 
postulated that only left-handed fermionic fields interact weakly, that is, the degree 
of parity violation is 100% ("V-A" interaction). This means that the SM corresponds 
to only g\ L being non zero, and leads to the following values of decay parameters 

p=l, 77 = 0, $ = 1, 8 = 1. (1.6) 

Many extensions of the Standard Model give rise to other couplings that mod- 
ify (HH). 

The contact interaction (jl.ip is not renormalizable, so only a tree level result can 
be obtained in a consistent way for the most general case. A calculation of radiative 
corrections requires either restricting the interaction to a V,A type, or specifying 
an underlying model leading to the effective interaction (jl.ip . Radiative corrections 
to the spectrum are significant |15j [T3"l [T6] and have been calculated under different 
assumptions by many authors. Very detailed results computed within the Standard 
Model are available 03 CE1 H [20l EE] . 



x + -5 | Ax 
3 



3 + 



Xn 



+ F^(x). 
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1.2 How general is the "general" matrix element? 

This section 1 discusses assumptions underlying (II. ip . 

First of all note that (II, lh assumes that all the four particles involved in muon 
decay are fermions. The two neutral particles, which are not observed in modern 
experiments measuring the muon decay spectrum, may be of a different nature. For 
example supersymmetric theories predict the decay of a muon into an electron and 
two light scalar sneutrinos mediated by a wino. Such a decay cannot be described 
accurately by a parameterization of the four-fermion interaction [7]. 

To understand what other assumptions are implied by (II. ip . let us consider the 
S-matrix element of the /i — > euu decay. There are four particles involved, fully 
described by 16 kinematic variables, the components of the four-momenta. These 
variables obey 4 relations pf = mf. The conservation laws corresponding to 10 
generators of the Poincare group impose 10 additional constraints. Thus 16—4—10 = 
2 independent invariants can be constructed from the four 4-momenta [22] . For 
example s = (p M — p Vp ) 2 and t = (p e — p u J 2 . In the general case the S'-matrix 
element should have the form [22] 

M fi =Y,fn(s,t)F n (1.7) 

I! 

where the functions of the kinematic invariants, f n , are called invariant amplitudes, 
and the F n are invariants which depend linearly on the wave amplitudes and 4- 
momenta of all the particles concerned. Equation (jl.7p may be understood to include 
all radiative corrections. 

The spin part of (jl.7p is more general than (jl.ip . The linear momentum de- 
pendence of F n in some cases can be absorbed into /„ by applying the equation of 
motion, for example: 

P2 (Hlal^iH^) = m 2 (^2}(^4) (1-8) 

But one can imagine a term for which this reduction will not work. For example, 
replace p2 — > P3 in (jl.8p . So, "derivative-free" is an assumption. This conclusion 
seems to contradict a statement in [23], p. 5, that "in the interaction term the 4x4 
differential operators can be reduced to (constant) 4x4 matrices." However, a con- 
tact interaction is assumed in [23] but not in the S'-matrix approach. So for the spin 

1 Material in this section is based on my term write-up for the Quantum Field Theory — II course 
at U of Alberta (2000). 
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part of the matrix element the assumption of "locality" of the effective interaction 
makes it also "derivative-free" . Extensions of the four-fermion interaction that allow 
for a linear momentum dependence of the spin part have been proposed \24\ [25] . 

The invariant amplitudes f n (s,t) in (jl.ip are just constants. To build f n we 
need a mass parameter, M, so that f n (s,t) = f n (s,t), where f n is a function of 
dimensionless variables s = s/M 2 , t = t/M 2 . Taylor expansion of the amplitude 
should look like 

/„W) = /„(0,0 ) + f ±- 2+ % £ + (1.9) 

In the context of a gauge theory M is the mass of an intermediate boson, e.g. M = 
Mw S> mn, and replacement of / n (s,£) by a constant can be justified. Again, this 
assumption is approximately equivalent to the assumption of a contact interaction: 
a heavy mediator means a short-range force. 

The interaction term (jl.ip contains also the assumption that lepton number is 
conserved. However, this assumption is not essential. Langacker and London [26] 
have shown that a Hamiltonian allowing both lepton flavor violation and total lepton 
number violation still leads to the same decay spectrum (|1.3jl . Moreover, there 
is a one-to-one correspondence between combinations of coupling constants of the 
lepton- number non-conserving Hamiltonian and coupling constants in (|l.ip . 

1.3 Tests for new physics 

with the muon decay parameters 

Early measurements of the muon decay spectrum helped to establish the current 
theory of the electroweak interaction. The much more precise experimental data 
available today are still in good agreement with the Standard Model. The best 
measurement of the muon decay parameter 5 before TWIST was [27] 5 = 0.7486 ± 
0.0026(stat) ± 0.0028(sys). In the rest of the chapter we will consider some con- 
straints on new physics that can be imposed by a more precise measurement of 
muon decay parameters. We will concentrate on parameters 5, a measurement of 
which constitutes the subject of this thesis, and £, which was constrained by the 
presented measurement of 5. At the end of this section we briefly mention some 
models where 5 differs from the SM value of 3/4, which have been excluded by other 
experiments. 
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1.3.1 Model-independent search for right-handed interactions 

By re-ordering (jl.2p we can write the fractions of decays where the muon interacts 
/i-handedly and the positron e-handedly as [7] 

QRR = \\9rr? + \9rr\ 2 (1-10) 

Qlr = \\g s LR \ 2 + \gU 2 + 3|<?Ll 2 (1-11) 
Qrl = \\9 S rl\ 2 + \9rl\ 2 + 3\9rl\ 2 (1-12) 
Qll = \\9 S ll\ 2 + \9 V ll\ 2 (1-13) 

The fraction of muons decaying through a right-handed interaction = Qrr+Qlr 
can be expressed through Michel parameters £ and 5 |14j : 

«-H 1+ i< -"«'}■ (L14> 

and thus is measurable by TWIST. The non-negative quantity Q R is exactly zero in 
the Standard Model. Any deviation from zero would indicate that the right-handed 
muon component participates in the decay process through either a scalar, or vector, 
or tensor, interaction. 

1.3.2 Left-right symmetric models 

In the Standard Model the charged weak current is purely V — A. A natural assump- 
tion is that the V + A current is suppressed, but not exactly zero |28J . Left-right 
symmetric models [29j [30j [311 [32j [33] extend the electroweak gauge group to include 
at least SU(2)r and refer to a spontaneous symmetry breaking mechanism to ex- 
plain parity violation. A general SU(2)l x SU(2)r x U(l) case is considered in [33] . 
The charged gauge boson fields are mixed: 

W L = cos(W 1 +sm(W 2 , (1.15) 
W R = (-sin(W 1 + cos CW 2 ), (1.16) 

where Wl, Wr are the interaction eigenstates, W±, W 2 , are the mass eigenstates, ( 
is a mixing angle, and to is a CP-violating phase. 

The Wr boson can contribute to muon decay only if the right-handed neutrinos 
are light enough, so that the process is kinematically allowed. The muon decay 
parameters affected in left-right symmetric models are p and £. Note that the 
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spectrum shape (jl.3p - (jl.5p depends on the combination P^£, not on £ itself. The 
polarization of muons from charged pion decay, P„, in left-right symmetric models 
is also different from unity. Introducing the notation 

2 2 

of, mi 

(L17) 

C, = — C (1.19) 
5L 



we can write 



P=^(1"2C 9 2 ), (1-20) 
C = l-2(t 2 + C 2 ), (1.21) 



»9 

P M = 1 - 2tj - 2C 2 - IteCg cos(a + w). (1.22) 



Here <?l, qr are the coupling constants, mi, rri2 are the masses of W\ and W2, 
V ud ' are the elements of the left- and right-handed quark mixing matrices, and 
a = &rg{V^ d } is a CP violating phase. (V^ d is chosen to be real.) 

It follows from (|1.2ip - (|1.22p that a measurement of constrains both the 
mass of the second charged gauge boson and the mixing angle. 

1.3.3 Non-local tensor interaction 

The ISTRA experiment observed a statistically significant deficit of tt~ —> e~v^ 
events in the P 7 > 21 MeV, E e > 70 MeV — 0.8P 7 region [35J. To explain it, a new 
momentum transfer dependent tensor interaction has been suggested |24| . This idea 
was discussed in the literature. In particular, [36] pointed out possible difficulties 
the hypothesis may have explaining nuclear beta decay data. However it could not 
be excluded [37] • Recently another experiment, PIBETA [38], also observed a deficit 
of radiative pion decay events (using 7r + ) in a similar kinematic region, renewing an 
interest in the problem. 

The suggested tensor interaction is non-local (momentum transfer dependent), 
and is not included in Eq. (II. ip . A new coupling constant, Qrr, needs to be intro- 
duced. The new interaction term, which should be added to (jl.ip . can be written 
as: 

-V2G F 9rr (eR\a a \\v e ) (^I^aI/Ar) (1-23) 
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This term contains only left-handed neutrinos and can interfere with the standard 
decay mode. This interference leads to a higher experimental sensitivity to that 
interaction. 

There is a field theoretical extension of the Standard Model [39] , which produces 
the effective tensor interaction fll.23|) . In [40] a prediction for the spectrum of 
positrons from muon decay is made based on the pion decay data. It is shown that 
the muon decay parameter 5 is very sensitive to the new interaction 



Standard Model value can be expected. 
1.3.4 Historical models 

The muon decay parameters have been discussed in the context of supersymmetric 
theories with light sneutrinos [3TJ [27] . However LEP data at the Z pole [32J 113 SU 
135] and above [46] I47| constrain rriy > 30 ... 94 GeV, depending on the assumptions. 
Therefore muon decay with sneutrinos in the final state is kinematically forbidden. 

An explanation of the LSND anomaly suggested by Babu and Pakvasa [38] in- 
volves a lepton number violation decay n + — > e + + VI + V\ (i = e, /i,r). Since the 
model requires p = d& 0.7485, it can be tested by TWIST. In 2003 the KARMEN 
collaboration put a strict limit on the emission of VI from fi + decay [39], excluding 
the explanation at 90% confidence level. 




(1.24) 



and with the suggested value g^ R 



0.013 almost a 10 3 deviation of 5 from the 
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Chapter 2 

Experimental setup 



The TWIST experiment is designed to measure the spectrum of positrons from muon 
decay in a wide range of energy and angle. A conceptual view of the spectrometer is 
shown on Fig. 12.11 An important feature of the detector is its planar geometry [50] . 
which gives the possibility to correct for the average energy loss of decay positrons 
with high precision in a data-driven way (chapter [6]) . The experiment uses a highly 
polarized surface muon beam |51j from the M13 secondary beam line at TRIUMF. 
The beam rate of about 2.5 x 10 3 muons per second is low enough to typically 
have no more than one muon at a time in the detector. A muon is stopped in the 
center of a symmetric stack of planar wire chambers and decays at rest. A 2 T 
uniform magnetic field preserves the direction of the spin of the stopped muons. 
The decay positron spirals in the magnetic field, leaving hits on the wires. The 
hits are recorded by TDCs and analyzed offline to reconstruct the trajectory of 
the particle and determine its energy and angle with respect to the magnetic field. 
A detailed description of the TWIST apparatus is given in The rest of this 

chapter summarizes different aspects of the experimental setup. 

2.1 Muon beam 

The TWIST detector is installed in the M13 secondary beam line [53] at the TRI- 
UMF cyclotron. Fig. 12.21 shows the M13 layout. The cyclotron produces a 500 MeV 
quasi-continuous proton beam, with 4 ns proton bunches striking a production tar- 
get every 43 ns. During the 2002 TWIST data taking a beryllium production target 
was used. Among the particles produced when beam protons interact with the target 
are positive pions. The dominant decay mode 7r + — > results in the production 

of muons, which can be transported through the M13 beam line to the experiment. 
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Beam pipe 

— Superconducting magnet and cryostat 




Figure 2.1: A drawing of the TWIST detector [52 




Figure 2.2: M13 beam line layout [53] . Bl and B2 are the dipole, and Q1-Q7 are 
the quadrupole magnets. The production target 1AT1 is seen by M13 at 135° with 
respect to the primary proton beam, the bends in Bl and B2 are 60° each. 
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The process tt + — > ^ + z/^ is a two body decay, thus the momentum of the /U + in 
the rest frame of the ir + is fixed, = 29.79 MeV/c. The relationship between the 
muon spin and the muon momentum is also predicted by theory. In the Standard 
Model the spin is antiparallel to the momentum. This relationship may be altered if 
the muon scatters in a material, since the Coulomb scattering of the nonrelativistic 
muons changes the momentum direction without influencing the spin. To preserve 
high polarization of the muon beam such interactions should be minimized. 

The surface muon beam technique [51] utilizes those pions that stop in the 
production target, then decay at rest. Passage through a material causes muons 
to also lose their momentum, and that loss, like the depolarization, is proportional 
to the amount of material crossed. By tuning the beamline to select muons that 
lost only a limited amount of momentum the depolarization can be controlled. The 
muons accepted by the properly tuned beamline come from pions decaying in a 
thin layer of material close to the surface of the target, with some "cloud" muon 
contamination from pion decays in flight. The 43 ns time structure of the beam 
makes possible the elimination of prompt particles produced at the time the protons 
hit the target, which includes the "cloud" muons. Since the life time of tt + is 26 ns, 
most of "surface" muons are emitted between the proton pulses. As discussed on 
page 1471 the measurement of S requires high muon polarization, though a knowledge 
of the precise value of the polarization is not important. For the 2002 TWIST physics 
data taking the M13 beam line was tuned to the momentum 29.6 MeV/c, with a 
momentum acceptance of 1.3%, resulting in a higher than 90% muon polarization 
seen by the TWIST spectrometer. 

At surface muon momenta the beam contains mostly positrons, muons, and a 
small fraction of pions [54] . The positron and pion beam backgrounds are removed 
by the reconstruction software (chapter [5]) . Data were also taken at 120 MeV/c for 
calibration purposes. At this momentum the beam predominantly contains pions. 

2.2 TWIST detector 
2.2.1 Wire chambers 

The TWIST apparatus uses wire chambers as the primary source of information. 
Two types of chambers are employed in the detector: drift chambers (DC), and 
proportional chambers (PC). They are similar in construction and use 15 /im sense 
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Figure 2.3: Side view of the TWIST cradle [52]. 

wires. Their cathodes are made from nominally 6.35 [im. thick doubly aluminized 
Mylar foils. (Monte-Carlo simulation of the detector, chapter U uses a measured 
mass density instead of the thickness, and also accounts for effect of stretching of 
the foils on their thickness.) The pitch of the sense wires is 4 mm for the DCs and 
2 mm for the PCs. The cathode-to-cathode distance is 4 mm in all cases. The 
PCs use a CF4/isobutane gas mixture, the high drift velocity of which provides fast 
response. One of the main functions of PCs in the analysis is to resolve tracks of 
different charged particles in time. The DCs use dimethylether (DME) gas, which 
gives high spatial resolution. 

The wires are positioned at 45° to the vertical direction to reduce the gravita- 
tional sag, therefore instead of the X and Y the chambers measure the U and V 
coordinates as defined in Appendix IIP. 41 The wire chambers are assembled into 
modules, each module having two or more wire planes. The volume between the 
chambers is filled with a helium(97%)/nitrogen(3%) mixture. A side view of the 
stack of wire chambers is shown on Fig. 12.31 The order, as seen by an incoming 
muon, of the modules is: PCs (4 planes), the "dense stack" (8 DCs), seven modules 
of "sparse stack" containing a pair of DCs each, the target module. The down- 
stream (after the stopping target) arrangement mirrors the upstream. The DC and 
PC planes are numbered sequentially, with the numbers increasing along the path of 
a muon. DC 22 is the last drift chamber and PC 6 is the last proportional chamber 
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Figure 2.4: Deviation of wires from their nominal positions, a = 3.3 /jm |52j . 

before the stopping target. 

The target module consists of 4 PC planes. The central cathode foil in this 
module also serves as the muon stopping target. During the 2002 data taking 
this target was made from 125 /xm Mylar, with conductive graphite coating on 
both sides. (A high purity aluminum stopping target was not available at the 
time.) Depolarizing interactions in the target rendered 2002 data unsuitable for an 
improved measurement of Pu£, but the extraction of 5 does not require a precise 
knowledge of the value of P M (page 1471) . 

An important advantage of the TWIST detector is the small amount of material 
in the tracking volume, leading to smaller effects from scattering and energy loss and 
thus to smaller related uncertainties. The thickness of one pair of DCs is only 1- 10~ 4 
radiation lengths. Also, the positrons cross only about 25 mg/cm 2 of material before 
entering the tracking volume at the first DC, compared to 240 mg/cm 2 before the 
tracking volume in the previous measurement [27] . 

A very high mechanical precision has been achieved in the production of the 
detector. The positions of the planes in the Z direction (along the beam) are defined 
by precise Sitall ceramic spacers with a negligible coefficient of thermal expansion. 
The cumulative error on the Z position is less than 50 /im over the whole length 
of the detector [52]. The positions of wires within each plane are accurate to a 
few microns, see Fig. 12.41 The relative alignment of different wire planes within 
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the detector has been accomplished using 120 MeV/c pion tracks taken with the 
spectrometer magnet off. The precision of this alignment is 5 fim (translations) and 
0.01 degrees (rotations) [52] . 

2.2.2 Electronics and DAQ 

A 195 //m thick scintillator mounted between the beam pipe and the upstream end 
of the TWIST detector is used to provide a trigger signal during physics data taking. 
The scintillator is read out by two PMTs and the trigger condition is a coincidence 
of the two signals. The thinness of the scintillator makes the system less sensitive 
to beam positrons, while muons, which have a higher density of energy loss at this 
momentum, produce a strong signal. An essential feature of the trigger is that it 
is unbiased: since the decay positron is not used by the trigger system, the trigger 
efficiency is not correlated with the muon decay parameters. 

An electric signal from a DC or PC wire is fed to a pre-amplifier mounted 
inside the wire chamber module. The output of the pre-amplifier is connected to 
a post-amplifier/discriminator in a CAMAC crate outside of the detector. The 
discriminator circuit provides a time over threshold signal, which is recorded by a 
multihit TDC with 0.5 ns time resolution. Upon receiving a trigger signal the TDC 
analyzes its internal buffer, and any activity from 6 fj,s before to 10 fis after the 
scintillator hit is read out via FASTBUS by a PowerPC. Time of the leading edge, 
as well as the width (time over threshold) of the signal are recorded for each of up 
to 8 hits per wire. Data are sent through an Ethernet connection to a dual 1GHz 
Pentium Linux computer running a MIDAS |55j based data acquisition system |56] , 
which writes them to a disk buffer, then to SDLT tapes. 

The gas gain of the drift chambers combined with the electronic amplification 
leads to an effective threshold of 1.6 electrons collected from a track to produce a hit 
|52j . The wire chambers operated at about 99.95% efficiency |52j . without a single 
dead or noisy channel during the 2002 data taking. 

In addition to the TDC data, the DAQ logs hundreds of "slow control" variables. 
They include voltages and currents for individual wire planes, gas flows through the 
chambers, proton beam current, temperatures at numerous locations in the TWIST 
detector, NMR measured magnetic fields of the spectrometer and beamline magnets, 
currents of the beam line magnets, atmospheric pressure, etc. 
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2.2.3 Spectrometer magnet 

The stack of wire chambers is placed inside of a 2 T superconducting solenoid. The 
solenoid, together with the outside steel yoke, produces highly uniform magnetic field 
in the tracking volume. The B z component of the field was mapped, Fig. 12.51 shows 
representative curves from the measurements. An OPERA-3d [58] simulation model 
was tuned to the measured B z . The OPERA-3d simulation produces the complete 
B(r) field (as opposed to simply B z (r)) that is used by TWIST Monte-Carlo and 
track reconstruction software. The simulated map reproduces the measured B z to 
better than 3 Gauss in the tracking region, giving the relative accuracy of 1.5 x 10 . 

2.2.4 Beam degraders 

To center the distribution of muon stopping position in the target, there is the 
possibility to fine tune the amount of material in the path of muons. This capability 
is provided by a gas degrader, a 21.67 cm long volume installed between the vacuum 
window of the beam pipe and the trigger scintillator. The gas degrader contains a 
He/C02 mixture. The fraction of CO2 can be varied from 0% to 100%, affecting 
energy loss of muons in the degrader and their final stopping position. 

It is also possible to install a plastic film in the path of muons. It was used to 
shift the stopping distribution to the upstream end of the stack of wire chambers to 
acquire the Monte-Carlo verification data (chapter [5]) . 
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Figure 2.5: B z vs z on the detector axis (top), and at the edge of the tracking volume 
(bottom). Note the zoomed vertical scale. The limits of the tracking volume in z, 
defined by the outermost DCs, are ±500 mm. The radius of the tracking volume is 
defined by the size of the wire planes. Plots from [57J. 
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Chapter 3 

TWIST data 



With a quasi-continuous beam from the TRIUMF cyclotron, the definition of a "data 
run" is arbitrary. In TWIST the DAQ was usually instructed to split data into files 
of about 1.9 GB each. One such file is a "data run". A typical run contains about 
8.5 x 10 5 data events (triggers), and was acquired in about 7 minutes for nominal 
surface muon beam. Data quality was monitored on per-run basis, and, if a problem 
was detected, a complete run was excluded from the analysis. 

A "data set" was defined as the amount of data required to achieve a statistical 
precision of ~ 10~ 3 on the muon decay parameters. While acquiring a data set, 
all controllable running conditions were kept unchanged. (However variations in 
e.g. atmospheric pressure could still introduce differences between runs within a 
data set.) Set A became significantly smaller than other sets because of an off-line 
rejection of bad runs. Table 13.11 summarizes data sets used for the extraction of 5 
and systematic studies. Note that set A has much lower statistics than other data 
sets. This is why set B was used in this work to quote typical numbers or show 
example plots. 
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Set name 


Dates 


Number 
of runs 


Fiducial 

events, 

millions 


Comment 


A 


Oct 8-9 


165 


7.9 


nominal 


B 


Nov 21-23 


318 


15.9 


nominal 


1.96 T 


Dec 2-4 


338 


16.5 


1.96 T spectrometer field 


2.04 T 


Dec 7-9 


240 


12.7 


2.04 T spectrometer field 


Cloud 


Nov 6-28 


561 


12.4 


Cloud muon beam 


DS Al 


Nov 29-30 


160 


7.7 


Outside materials systematic 


Slightly 
Upstream 


Oct 5-8 


307 


7.3 


Stopping location systematic 


Low rate 


Oct 13-20 


338 


17.9 


Beam intensity systematic 


High rate 


Oct 11-13 


341 


14.1 


Beam intensity systematic 


B2 + 10G 


Oct 20-21 


348 


15.4 


Channel magnets systematic 



Table 3.1: Data sets mentioned in the thesis. Fiducial region is defined in chapter [5j 
All dates are in 2002. 
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Chapter 4 

Monte-Carlo simulation 



TWIST uses a very detailed Monte-Carlo simulation program (MC), which is based 
on the GEANT [59] package. It produced digitized output in the same format as 
the DAQ, except that additional MC-specific information may be included. Pro- 
duction of a large amount of Monte-Carlo events, matching TWIST data statistics 
(chapter [3]) , was made possible by the use of WestGrid computing facility |60j . 

The geometrical description of the detector contains all the components of the 
hardware with which a muon or a decay positron could possibly interact. Each 
individual wire of the wire chambers is implemented in the software. The wire 
planes are offset and rotated to their as-measured positions. A map of the magnetic 
field (sections I2.2.3|) . which extends to the outside of the yoke, is used to propagate 
charged particles in the simulation. 

The initial kinematics of an event contains a muon, and possibly other muons 
and/or beam positrons. The probability of having the pile-up particles is determined 
by the specified muon and beam positron rates. The positions and directions of flight 
of the muons are sampled from experimentally measured distributions {x,dx/dz} 
and {y,dy/dz} and reproduce the observed position-angle correlations. The beam 
particles, muons and positrons, are started outside of the yoke, and GEANT track- 
ing propagates them through the fringe field of the spectrometer magnet into the 
TWIST detector. 

In the TWIST Monte-Carlo, unlike the standard GEANT3, the direction of the 
muon spin is also tracked in the magnetic field. The initial spin direction is defined as 
antiparallel to the muon momentum. Depolarizing interactions of the stopped muons 
are simulated as a step function followed by an exponential relaxation. Setting the 
initial polarization to —0.935 and the time constant to 5.8- 10 -5 s [61] reproduces the 
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behavior of P^(t) in data as it is observed for t > 1 fj,s (this time limit is discussed 
in chapter [5]). A mismatch in polarization between data and Monte-Carlo would 
not bias the value of 5 in the analysis, but a large difference in the average 
could require generating larger "derivative" (chapter [7|) samples to achieve the same 
statistical precision. 

The energy loss of a muon is simulated by GEANT, and determines the stopping 
position of the particle. The simulation of the muon stopping process has been val- 
idated using special data runs with muons stopping in the middle of the upstream 
part of the chamber stack. Because each of the chamber modules is much thin- 
ner than the stopping target, the muon stopping distribution in these runs spreads 
out over several wire chambers, and the distribution of the last (the most down- 
stream) wire plane hit by muon can be used to observe the shape of the stopping 
distribution. It has been shown [62] that the simulation matches the shape of the 
stopping distribution well, but a constant offset equivalent to an about 86 //m of 
additional plastic (Mylar) is required in Monte-Carlo to match the mean stopping 
position of the muons. The peak of the stopping distribution within the target is 
not directly observable in data, but the tails of the distribution are still accessible 
through the last plane hit information. The following procedure was used to de- 
termine the setting of the gas degrader for the nominal data taking. A histogram 
of the last muon hit from Monte-Carlo, with the stopping distribution centered in 
the target, was compared with similar histograms from data for different settings 
of the gas degrader. The setting corresponding to the best match to the Monte- 
Carlo distribution was used for the data taking. All nominal Monte-Carlo sets were 
generated with the muon stopping distribution centered in the target. The energy 
calibration procedure (chapter [6]) compensates for any remaining differences in the 
average muon stopping position. 

A muon decay subroutine returns the energy and angle of the decay positron with 
respect to the muon spin as dictated by the Michel parameters input. The theoretical 
decay spectrum includes full O(a) radiative corrections with exact electron mass 
dependence, as well as leading and next-to- leading logarithmic terms of 0(a 2 ), 
leading logarithmic terms of 0(a 3 ), corrections for soft pairs, virtual pairs, and 
an ad- hoc exponentiation [17\ \TE[ fT9j [20} 121] . The actual implementation of the 
muon decay spectrum is separated from the rest of the Monte-Carlo code to make 
possible a blind analysis, as is explained in section 17.41 
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The wire chamber response is simulated by randomly creating ionization clusters 
along the path of a charged particle when it crosses a drift cell, calculating the drift 
time of each cluster to the wire, and simulating the overlap of different clusters to 
produce an above-the-threshold signal. The obtained time of the signal is smeared 
to simulate electronics effects. The parameters of the method are derived from a 
detailed GARFIELD [63] simulation study [63], and the "electronics" smearing is 
tuned to TWIST data. The critical piece of information affecting the accuracy of 
track reconstruction (chapter [5]) is the time of the leading edge of DC chamber 
signals. A good match of the simulated distribution to TWIST data has been 
demonstrated |65| . 

Full GEANT physics interactions are enabled, so such processes as creation of 
delta-electrons, or conversion of a bremsstrahlung photon into an electron-positron 
pair, may create additional hits on the detector wires. 

Interactions of decay positrons in detector materials distort the reconstructed 
spectrum and can lead to biases in the values of the measured decay parameters. 
TWIST relies on Monte-Carlo to compensate for these effects (chapter [7]) . Know- 
ing the accuracy of the simulation of the interactions is important to determine 
the corresponding systematic uncertainty (chapter [8]). To verify a claim that in 
GEANT3 

the cross-sections of the electromagnetic processes are well reproduced 
(within a few percent) from 10 keV up to 100 GeV, both for light (low 
Z) and for heavy materials 

(|59|. section PHYS001), data runs with muons stopped in the upstream end of the 
detector have been taken. In these runs a decay positron emitted in the downstream 
direction crosses the whole spectrometer. The upstream and downstream halves of 
the detector can be used as two independent devices to measure the momentum and 
angle of the positron before and after the central target, thus allowing an extraction 
of the energy loss and angular scattering distributions from data. These distributions 
can be compared with similar distributions obtained from a corresponding Monte- 
Carlo simulation. The Monte-Carlo validation data were taken in 2003, with a high 
purity aluminum stopping target of a known thickness. Thus the validation, unlike 
the final result, is not affected by the uncertainty in the thickness of the graphite 
coating on the Mylar stopping target used for the main 2002 data sets. 
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Figure 4.1: Positron momentum loss in data (solid) and Monte-Carlo (empty mark- 
ers). Top: linear vertical scale, bottom: logarithmic scale. The Monte-Carlo his- 
togram is normalized to data. Mean values of the distributions are —126.9 keV/c 
(data) and -122.9 keV/c (MC), the difference is 4.1 ± 1.3 keV/c. The RMS is 
0.269 keV for data, 0.258 keV for MC. (The RMS in this study does not represent 
TWIST momentum resolution, see text.) Analysis by Rob MacDonald. 
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Figure 4.2: Positron angle change in data (solid) and Monte-Carlo (empty mark- 
ers). The Monte-Carlo histogram is normalized to data. The mean angle change is 
—2.9 mrad for data, —1.9 mrad for MC. The RMS is 16.4 mrad for data, 16.6 mrad 
for the Monte-Carlo. (The RMS in this study does not represent TWIST angular 
resolution, see text.) Analysis by Rob MacDonald. 
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Fig. l4.1l shows an overlay of the momentum loss distributions for data and Monte- 
Carlo. There Ap = ^downstream — Pupstream is the difference between the positron 
momenta reconstructed by the two halves of the detector. In TWIST geometry mo- 
mentum loss of a positron is proportional to 1/| cos(0)| (see chapter [6|); a factor of 
| cos(#)| is used in the study to compensate for a possible difference in the angular 
distributions of the emitted positrons between data and Monte-Carlo while com- 
paring the momentum losses. The difference between the data and the simulation 
does not exceed 5% for both the average losses, and the widths of the distributions. 
To separately look at "hard" processes, an arbitrary bound of 1 MeV/c was de- 
fined. The tails of the distributions shown on Fig. 14.11 were integrated from — oo to 
— 1.12 MeV/c (i.e. to 1 MeV/c below the average). The discrepancy in the fraction 
of such "hard scatter" events was found to be 14%. Therefore the numbers we used 
to estimate a systematic uncertainty due to the quality of GEANT simulation of 
positron interactions are 14% for hard interactions, and 5% for intermediate and 
soft interactions. 

The distribution of AO = ^downstream — ^upstream is shown on Fig. I4.2L The width 
of the distributions agrees to about 1%, but there is a noticeable shift in the mean 
value of AO. The interpretation of the mean value of this distribution is complicated. 
Since there is more phase space at higher angles due to the sin 6 factor, a naive 
expectation is to observe a small positive shift in A9. This has been confirmed by 
a Monte-Carlo study [66j, which demonstrated a ~ 1 mrad positive bias. An effect 
of a non- uniformity of the magnetic field is negligible [66] . However the distribution 
of the reconstructed AO has a negative mean for both data and Monte-Carlo. Thus 
is should be attributed to biases in the track reconstruction, which affect the two 
halves of a track in a different way. At least two causes for such biases are known. 
In the GEANT validation studies positrons originate in PC 4, and no information 
from PC 1-3 is available for pattern recognition (chapter [5]) for the upstream part 
of a track, but full information is available for the downstream part. Also, the track 
fitting always assumes that a positron originates in the stopping target to determine 
the sign of a time-of-flight correction to drift chamber hits. Therefore the correction 
was applied with the wrong sign for the upstream parts, and with the correct sign 
for the downstream parts of tracks in the validation study. Such biases are common 
for the data and the simulation, and may shift the distributions to negative AO. A 
possible explanation for the difference in the central values between data Monte- 



24 



Carlo is a misalignment of the detector to the magnetic field. Chapter [8] evaluates 
a systematic uncertainty associated with such a misalignment. 

It has to be noticed, that the widths of the distributions on Fig. 14. II and Fig. 14.21 
do not represent TWIST resolutions. The lack of information from PC1-3, and 
the wrong time-of-flight correction mentioned above worsen the tracking quality. 
Moreover, a positron in these studies goes through twice as much material as a 
positron from a muon stopped in the target. Also, a track is reconstructed in the 
two halves of the detector, so random fitting errors contribute twice. 
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Chapter 5 

Spectrum reconstruction 



The primary purpose of the spectrum reconstruction is to produce a 2-dimensional 
spectrum of decay positrons in momentum and angle using information available 
from the TWIST detector. 

Real and simulated data are analyzed by the same reconstruction code in es- 
sentially the same way. The only non-trivial difference in treating the two cases is 
the crosstalk removal, which is discussed below. Technically the analysis was done 
in two stages. The first stage, required a large amount of computations and was 
performed at West Grid [60]. It consisted of the following steps: 

Crosstalk removal. A hit on a wire may induce a signal in a different channel 
through electronics crosstalk. Most of crosstalk in TWIST is induced by highly 
ionizing muons, therefore it does not affect decay positrons a microsecond later. 

The effect has been studied in detail in the hardware by pulsing a channel 
and observing the crosstalk signal on an oscilloscope. Crosstalk signals have 
characteristics allowing for their identification and removal by software |67j : 
there is a 5-65 ns delay from the generating pulse, and the width of a crosstalk 
signal is much smaller than the width of a real signal. 

An algorithm utilizing these characteristics removes crosstalk hits from real 
data events before any other analysis is done. Because crosstalk is not simu- 
lated by TWIST Monte-Carlo, and also because the width of chamber signals 
in the simulation is not tuned to data, this algorithm is not invoked on Monte- 
Carlo events. 

Windowing. Chamber hits are grouped in time. Different groups correspond to 
tracks from different particles. 
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Classification. For each time window, characteristics of its chamber hits are used 
to guess what process has occurred there [68] . For example, a decay positron 
hits wire chambers in only one half of the detector, while a beam positron 
crosses the whole detector. Based on the number and type of time windows in 
an event, and the times of different windows, a classification code is assigned 
to the event [68] . 

Pattern recognition. A cluster of consecutive wires hit in a plane determines only 
one transverse coordinate of a passing track. Information from an orthogonal 
pair of wire planes can be used to determine a "space point" on the track. 

For a time window, the pattern recognition uses space points formed by the 
hits belonging to the window to find parameters of possible helical tracks 
passing through the points. It uses a combinatorial technique [69] 170] and can 
find more than one track per window. 

A feature of the arrangement of the wire chamber modules in the detector, 
which can be seen on Fig. 12. 3\ is that the distances between the corresponding 
planes of different pairs form a pattern: 5.2 cm, 7.2 cm, 5.2 cm, 7.2 cm, 5.2 cm, 
7.2 cm, with only two independent distances, for the 7 modules closest to the 
stopping target on either side. That means that a helix with the wavelength 
L = 5.2 + 7.2 = 12.4 cm has only two distinct measurements of its transverse 
position, so its radius can not be reconstructed using the space points [7l~|[72]. 
(In a projection along the detector axis, the 7 space points collapse into only 
2 points on the circle — projection of the helix.) 

The "dense stacks" of wire planes at the outer edges of the detector, along 
with information from PC chambers, help to resolve the ambiguity. However 
the worsening of the quality of reconstruction of such tracks still prompted 
the introduction of a fiducial cut on p z , described below, to stay away from 
the "magic wavelength" zone. 

Wire center fits. Tracks are fit to the positions of hit wires in a time window. The 
"narrow windows" technique [73] is used. To account for multiple scattering, 
kinks at the positions of sparse stack chamber modules are allowed [71] . The 
magnetic field map (section I2.2.3P is used. The resulting fit has sufficient 
precision to resolve most left-right ambiguities. 
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Track fitting The parameters of the tracks are refined using drift time information. 
The necessary space-time relations are obtained from GARFIELD. Kinks and 
magnetic field map are used. 

No cuts on events were imposed at the first stage, and a summary of each event 
was written out in a "ROOT tree" format [7S]. The output data were subsequently 
analyzed to select events and tracks to be used in the final spectrum. The following 
cuts, illustrated on Fig. 15.11 were used. 

TCAP. This is a cut on the time of flight of the trigger particle through the M13 
beam line. It selects muons from stopped pions that decayed in between the 
proton pulses from the cyclotron, and rejects "cloud" muons (section l2.ip . The 
purpose of the cut is to improve muon polarization. It also rejects triggers from 
beam pions and prompt positrons. 

The time of flight cut is only applied to data events. The Monte-Carlo does 
not simulate the 43 ns time structure of the beam, and is instead produced 
using measured parameters of the surface muon beam with the TCAP cut 
applied. 

Event Type. This cut uses the event classification code, and is a combination of 
several requirements. The event must be triggered by a muon, and not a 
beam positron. There must be only one identified muon. A unique decay time 
window must be identified. 

The decay must happen at least 1.05 /is after the trigger to ensure that DC hits 
with the longest drift time ~ 1 /is from the muon do not affect reconstruction 
of the decay positron. Since the muon life time is about 2.2 fis, this requirement 
alone rejects 38% of events. 

Pile-up beam positrons are allowed, but they must be well separated in time 
(at least 1.05 //s) from both the muon and the decay positron to avoid an 
overlap of DC hits. Events with DC overlap constitute another major fraction 
of all events rejected by the Event Type cut. It is important to note that 
rejection of DC overlap events does not introduce a bias in the spectrum, 
because the probability of an overlap does not depend on the momentum and 
angle of the decay positron. 

On the other hand, beam particles within ~ 100 ns of the decay can not be 
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Figure 5.1: Top: each bin contains the number of events before the cut; bottom: 
the number of events rejected by cut. Solid markers: data, empty markers: MC. 
Monte-Carlo is normalized to data by matching the number of accepted events. 
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reliably separated by the classification code from such processes as creation 
of delta-electrons or backscattering of the decay positron. Therefore all such 
"PC overlap" events are kept at this stage. 

Muon Z. Require that the event is consistent with the muon stopped in the central 
target, i.e. require that the last chamber hit by the muon is PC 6. 

Muon radius. Constrain muon stopping position in the transverse plane to r < 
2.5 cm. 

Decay Window Time. Require the decay to occur before 9 /xs from the trigger, 
so that long drift time hits are not cut off by the 10 fxs DAQ limit. 

N tracks Require at least one track candidate, as defined by the pattern recogni- 
tion, in the decay window. 

All the cuts above are essentially event-level cuts. An event may contain more 
than one successfully reconstructed track. There are several reasons for multiple 
tracks to be found. For example, due to a hard scatter a decay track may be split 
into two helical parts, before and after the scatter, by the pattern recognition. That 
would lead to two fit tracks for a single decay positron. Or a decay positron may 
backscatter off material outside, re-enter the tracking volume, and produce a second 
track, which may cross the whole detector. Such events may be indistinguishable 
from a decay being overlapped by a beam positron particle. A genuine beam positron 
overlap is yet another possibility to produce multiple tracks. 

The challenge is to identify the "correct" track to be included in the spectrum. It 
is handled as following [76] . A set of "decay candidate tracks" is created. Initially it 
consists of all tracks identified in the decay window by the pattern recognition. Then 
each of the candidates is subject to the cuts below. Failed candidates are eliminated 
from the set. A corresponding event-level cut is defined as the requirement that the 
set of decay candidates is not empty after the track cuts. 

ierror Eliminate track candidates that did not produce a successful fit. 

startstop Require that the ends of a fit track are on the "correct side" of the 
detector, as determined by the classification. I.e. the track begins and ends 
in the upstream half for upstream decay type events, in the downstream half 
for downstream decay type events. Tracks that cross the central target are 
always rejected here. 
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charge The direction of the winding of the helix must be consistent with a positive 
particle that originated in the central target. 

pair matches This piece of code essentially attempts to "glue" pieces of tracks 
that were split in the reconstruction because of a hard scatter |76j . 

The cut is done in two stages. First, for each track among the decay candidates, 
a set of "anti-tracks" is found as described below. An "anti-track" for a decay 
candidate is a track in the decay window that, if successfully "glued" to the 
candidate, would indicate that the candidate track should not be used. 

Then the closest distance of approach 1 (CDA) is calculated between the can- 
didate and each of its anti-tracks. If CDA is less than 0.5 cm for at least one 
anti-track, the candidate is rejected: probably the anti-track and the candidate 
track belong to the same particle. The decision is done on purely geometrical 
grounds: do the two tracks intersect? They can belong to the same particle 
only if they do. 

Selection of anti-tracks 

The set of anti-tracks for a decay candidate is found through the following 
procedure: 

1. Start with all good fits in the decay window, but exclude the decay can- 
didate itself, to form a set of anti-track candidates. 

2. Exclude all tracks which overlap with the decay candidate in Z. If two 
tracks have hits in the same DC plane, they can't be from the same 
particle. (Sharing of DC hits in not allowed in the track fitter. One DC 
hit could belong to no more than one fit track.) 

3. Exclude all tracks which are in the same half of the detector and farther 
from the target than the decay candidate. We want to keep the decay 
candidate because it is closer to the target, and therefore provides more 
accurate information on the momentum and angle of the positron at the 
decay point, even if the anti-track belongs to the same particle. 

At this point an event is accepted and will produce an entry in the final spectrum. 

However there is still no guarantee that a unique decay track has been identified. 

1 The distance used is the shortest distance between two tracks in the transverse plane, not in 
3D. 
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Figure 5.2: The mean number of extra decay track candidates before cut or select. 
Solid markers: data, empty markers: MC. 

In fact about 0.7% of events still contain more than one decay track candidate, as 
can be seen on Fig. 15.21 The following algorithm [76] is used to decide among the 
multiple candidates. These "selects", unlike cuts, never reject all candidates, but 
choose one. 

dplane to target Find the DC plane closest to the central target, which is used 
by a decay candidate. Eliminate all candidates that do not use a hit from that 
plane, and therefore start farther from the target. 

mu-e vertex The remaining < 1CP 4 fraction of multiple track cases is resolved 
by the proximity of extrapolation of the positron tracks to the muon stopping 
position. Since the error on the two muon coordinates is significantly different 
because of the large scattering near its stopping position, an elliptical metric is 
used: i? c iii P sc = (U e - U^) 2 + (V e - V^/a^, with a vu = 1.7 found empirically 
by comparing RMS of muon-positron mismatch in the U and V directions. 

The code is guaranteed to select only one decay track per event. 

The momentum and cos(#) values of the selected tracks represent an unbinned 
decay spectrum. This "raw" spectrum was used to perform the energy calibration 
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procedure, described in chapter [6j Then calibration results were applied to the 
unbinned decay spectrum (chapter [6]), and the corrected spectrum was filled into 
a 2-dimensional histogram in momentum and cos(#), which was further used to 
perform a fit and extract values of the muon decay parameters (chapter [7J. An 
example of a reconstructed data spectrum is shown in Fig. 15.31 

The fiducial region used in the extraction of 5 is defined by the following con- 
straints: 

p < 50 MeV/c The shape of the reconstructed spectrum at the end point is de- 
fined by the detector resolution. On the other hand, in the bulk of a smooth 
spectrum its distortion due to detector resolution is a second order effect. Ex- 
cluding the end point region from the decay parameter fits (chapter [9|) drasti- 
cally reduces the sensitivity of the result to discrepancies between resolutions 
for real and simulated data. 

Another reason to stay away from the end point is to keep the decay parameter 
fits statistically independent from the energy calibration fits (chapter [6]) . 

\p z \ > 13.7 MeV/c This cut eliminates tracks affected by the "magic wavelength" 
problem, which is discussed above in this chapter. 

Pt < 38.5 MeV/c That requirement, together with the r < 2.5 cm cut on the 
muon stopping position, insures that the decay positron track is radially con- 
tained within the instrumented region of the detector. 

0.50 < | cos#| < 0.84 Events at high angles (small | cos(#)|) are more affected by 
multiple scattering and momentum struggling, leading to worse resolution. 
They are also more difficult to reconstruct. 

At small angles, it becomes difficult to determine the wavelength (p z ) of the 
tracks. Reconstruction biases observed in this region lead to a deviation of 
the average energy loss p rcc — pmc from the 1/| cos(#)| behavior required for 
energy calibration (chapter [6|) . 

The shape of the fiducial region in the momentum and cos(#) plane can be seen 
in the upper left panel of Figs. I9.1H9.5I 

The average resolutions in the fiducial region are: 150 keV (FWHM) for momen- 
tum, 0.015 rad (FWHM) for the 6 angle, and 0.01 (FWHM) for cos(0). Distributions 
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of the differences of reconstructed and Monte-Carlo values for events reconstructed 
within the fiducial are shown on Fig. 15.41 
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Figure 5.3: Reconstructed data spectrum from set B. a) In the momentum-cos^) 
plane, b) Momentum spectrum for 0.50 < |cos(0)| < 0.84, proportional to F\s(p). 
c) Difference of momentum spectra for —0.84 < cos(#) < —0.50 and 0.50 < cos(#) < 
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Figure 5.4: Distributions of differences of reconstructed and Monte-Carlo values. 
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Chapter 6 

Energy calibration 



The momentum of a decay positron measured by TWIST is affected by energy loss in 
the detector materials. For a muon stopping distribution that is not centered in the 
target the average momentum shift is different for upstream and downstream decays, 
introducing an asymmetry in the spectrum. Also the average energy loss may be 
different in data and the simulation even for a centered stopping distribution. The 
momentum also scales with the ratio of the "true" magnetic field seen by the particle 
to the field used by the reconstruction program. The ratio is unity for the Monte- 
Carlo, but the field map is not a perfect representation of the real field. These effects 
may lead to different spectrum distortions for real and simulated data, producing a 
bias in the extracted values of the muon decay parameters. The goal of the energy 
calibration procedure is to correct the reconstructed spectra to compensate for these 
differences between the data and the simulation. 

The calibration is done on physics data, using the same reconstructed spectrum 
as for the extraction of the decay parameters. The calibration point is provided 
by the sharp edge of the muon decay spectrum at the upper kinematic limit. Its 
position is determined by the muon and the positron masses, and is therefore known. 

The planar geometry of the TWIST spectrometer leads to an exact 1/| cos(#)| 
dependence of the amount of material traversed by a decay positron [50]. This 
dependence is obvious for the flat stopping target, cathode foils, and the gas layers. 
For the cylindrical wires, that dependence still applies on average, because the 
probability of hitting a wire changes as l/|cos(0)|. In the 15-53 MeV/c TWIST 
range of momenta, this translates into a linear dependence of the average energy loss 
of the positron on 1/| cos(#) | . That known dependence provides a way to disentangle 
the effects of the magnetic field scale, which is angle-independent, and the energy 
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loss. The position of the spectrum end point as a function of angle is given by 

p * w = ( 1+ ^)( a -S»l)' (<u) 

Here p r e ^ e (d) is the end point position of the "raw" reconstructed spectrum, po ~ 
52.828 MeV is the kinematic spectrum limit, and the constants a and (3 describe the 
momentum loss and the field scale mismatch, respectively. To accommodate a non- 
centered muon stopping distribution, the calibration procedure allows for different 
energy loss parameters in the upstream and the downstream parts of the detector, 
a = a u or a = a^- The difference of the upstream and downstream energy losses 
oyiff = a u — c*d is proportional to an offset of the muon stopping distribution from 
the centered position, while the sum a sum = a u + does not depend on the muon 
stopping position and is a measure of an effective thickness of the detector. 

A determination of the calibration parameters /3, a u , ad is not trivial. The chal- 
lenge is to determine the end point position, p T c ^ c (0), on reconstructed data, where 
the sharp edge of the theoretical spectrum is smeared by the detector resolution. 
A model function describing the spectrum shape at the end point is needed. This 
complicated shape is a convolution of the muon decay spectrum and the detector 
resolution, and a model function can not be expected to perfectly describe the data 
distribution. Therefore the result of a fit depends on the range of momenta used. 
An objective procedure establishing the range must be developed. Since the bias 
of a fit depends on the relative position of the fit range with respect to the end 
point, that procedure should be adaptive and produce the same relative position 
for different absolute positions of the end point in the input data. It also has to 
be noticed that fitting in the end point region by definition involves transition from 
high to low per bin statistics, and care must be exercised in handling the statistical 
issues properly. 

A straightforward approach of fitting p r e ^ e (0) independently for different angles, 
with an adaptive choice of fit range, and then applying Eq. (|6.1|) to fit a straight 
line through the resulting points, has been tried \77\ I78j. Different end point model 
functions were used, some including the effects of radiative corrections on the muon 
decay spectrum [77]. An important result of these studies is a demonstration that the 
end point indeed has a 1/| cos(0)| dependence. However to stabilize the independent 
end point fits, especially for small downstream angles where the statistics is low, a 
large range of momenta ~ 2 MeV/c needs to be used. It is more difficult to find 
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Figure 6.1: A convolution of the cut-off linear function with a Gaussian gives the 
shape on the right, which is used to fit the end point of the reconstructed spectrum. 

a model function describing the data in a large range. Another issue is that the 
shape of the data distribution becomes sensitive to the values of the muon decay 
parameters if a large range of momenta is looked at, which is an undesirable effect. 
Also, there is no single goodness of fit criteria in this approach. The goodness of 
the final straight line fit does not include any information about goodness of the 
individual end point fits. 

To overcome these issues, a fitting procedure that uses a global fit to the 2-dimen- 
sional momentum and cos((9) reconstructed distribution has been developed |79j . 
The model function is constructed as following. The momentum dependence is 
given by a convolution of a "slope and a step" shape 

a(l + by)Q(y), where y = p - p edge (6.2) 

with a Gaussian, illustrated on Fig. I6.ll Here Q(y) is the Heaviside step function. 
There are 4 parameters a, b, p c d gc; which defines the position of the end point, and 
the Gaussian a. The explicit form of the end point model is 

\{a + by) erfc(^) - -|= «p(-^). (6.3) 

This one-dimensional function is used to describe the momentum dependence of the 
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spectrum for a fixed angle 9. To obtain a 2-dimensional function, a dependence of 
the parameters a, b, p e dge and a on 9 is introduced. For p e dge 5 the expected linear in 
1/| cos(0)| behavior from Eq. (16. If) is used, with different momentum loss parameters 
a u and but common field scale /3 in the upstream and the downstream: 

raw ( ^ _ (, , P \ ( n ^e(-COs(g)) + q d e(cOs(g)) \ 
Pedge W - ^1 + -) [PO ) ■ (6.4) 

The muon decay spectrum is linear in cos(0), therefore a linear parameterization 

a(9) = ao + a\ cos(#) (6-5) 

is used for the normalization. For a and b, suitable parametrizations were found 
empirically by fitting Eq. (16. 31) to the data distribution independently for different 
angular bins, and observing the behavior of the fitted parameters as functions of 
angle. That resulted into the following choice: 

6(0) = 6 o + 6 lC os(0). (6.6) 
°W = |-^, (6.7) 



Equations (|6.3p ()6.7h completely define the shape of the fitting function. They 
contain 8 parameters: j3, a u , ad, ao, a\, bo, b\, ctq, that are determined from a fit 
to a reconstructed spectrum. Fig. 16.21 shows an energy calibration fit to data for 
several angles. 

The fit is done by minimizing —2 In A, where A is the Poisson likelihood ratio [80] : 

N 



- 2 In A = 2 ^2 



i=l 



f 1 Ui 

fi-rn + iii In — 

Ji . 



(6.8) 



In (|6.8p the summation runs over the bins of the histogram in the fit range, rtj is the 
number of entries in a bin, and /j is the expectation value for rii computed using 
Eqs. (|6.3p - (|6.7p . The advantage of the binned likelihood statistics —2 In A is that it 
can be used not only to perform the fitting, but also to easily estimate a goodness 
of the resulting fit [80] . In the large sample limit the minimum value of —2 In A is 
distributed as \ 2 - 

As is mentioned above, for a fixed angle the same momentum fit range relative 
to the end point must be selected, so that no biases between data and Monte-Carlo 
are introduced due to different absolute positions of end points. A series of tests, 
described in Appendix I10.4| led to the following choice: 

Ped g e(#) - 0.75 MeV/c < p < Pcdgc (9) + 0.5a(9). (6.9) 
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The fit is done iteratively. It uses its current estimate of p e dge(#) to fix the fit range 
for the next iteration. There is no explicit limit on the number of iterations. Instead 
each new 2-dimensional fit region £1/. is compared with all previously seen regions 
Qi, . . . , fifk—l- The stopping condition is £1^ = ^k-i- Because a binned fit is used, a 
rounding to histogram bin boundaries ensures that the number of possible fit regions 
is finite, so the fit always converges. In practice it usually takes < 10 iterations. 
It is common to have I > 1, in that case the average value of /3 over iterations 
k — l + 1, . . . ,k is computed, and the iteration with /3 closest to the average is chosen 
as the final result. This is an arbitrary procedure, but the spread of the results over 
the regions of the "convergence cycle" is always much smaller than the statistical 
error of the fit. 

Table 16.11 shows results of energy calibration fits for data sets and correspond- 
ing Monte-Carlo sets used in the extraction of 5. A systematic difference in the 
momentum resolution parameter gq at the end point is apparent from the numbers. 
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Table 6.1: Energy calibration results. All parameters are in keV/c. The large 
deviations in f3 for two of the Monte-Carlos are because of a mistake in setting the 
field scale in MC production. The energy calibration procedure corrects for this 
mistake. 

This discrepancy may arise for multiple reasons. For example, the spectrum re- 
construction uses the same drift chamber space-time relations for data and Monte- 
Carlo. These relations may be slightly different in the real detector, but the simula- 
tion uses exactly the same drift tables as the reconstruction to generate Monte-Carlo 
events. Therefore the reconstruction of Monte-Carlo can be expected to perform 
better. This holds true for all other calibrations, such as alignment corrections, 
electronic timing offsets TO, etc. The GEANT handling of positron interactions in 
detector materials may also contribute to the effect. All the calibrations are counted 
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Figure 6.2: The end point of the muon decay spectrum, and sections of the 2- 
dimensional end point fit function for several angular slices. The smallest and the 
largest angles in the upstream (top), and the downstream (bottom) are shown. Data 
set B. 
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as sources of their corresponding systematic uncertainties, as well as an imperfection 
of the GEANT simulation, chapter [HI However an effect of the momentum resolu- 
tion discrepancy on the final result, regardless of its cause, has been estimated and 
included in the final estimate of systematic uncertainty, chapter [8j 

Typical correlation coefficients among fit parameters are shown in Table 16,21 
Since the parameters of interest /3, a u and ad are highly correlated, their correlation 
must be taken into account while estimating a systematic uncertainty of the result 
due to energy calibration. 
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Table 6.2: Correlation coefficients for energy calibration fit to set B. 



An energy calibrated decay spectrum is obtained from the corresponding "raw" 
spectrum by recomputing the reconstructed momentum for every event: 



P 



P 



+ 



a 



(6.10) 



1 + P/po | cos (6>) | 

The momentum is scaled by the magnetic field scale, and corrected for the energy 
loss, a = a u for upstream decays, ay for downstream, po is defined after Eq. (|6.ip . 
This calibration procedure brings the end point of the calibrated spectrum to its 
"theoretical" value, p e dge(G) — Po- Checks by running energy calibration fits on 
calibrated (instead of raw) spectra were done. As expected resulting (3, a u , and ad 
were consistent with zero. 
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Chapter 7 



Method for extraction of the 
decay parameters 



A reconstructed spectrum differs from a theoretical one because of finite experi- 
mental resolution and because the interaction of the decay positrons with detector 
materials changes the energy and angle of the particles. 

Consider a parametrized theoretical probability distribution function /theor(^'; A), 
x' € f^o, where A are the parameters of the theory, and the x' are the kinematic 
variables of interest. The probability K{x, x') dx of reconstructing an event that 
occurred at x' in a volume dx around some point x defines the response function 
K, which describes the combined effect of the detector and the reconstruction. The 
reconstructed spectrum / rec (x;A) can be written as 



where b{x) is the background term. 

Several approaches can be used to deduce the theory parameters A from a mea- 
sured spectrum, / rec . One can try to solve (|7.ip for /theor- Deconvolution methods 
are available [81], and were used by some experiments (e.g. [82]). However they are 
not practical for TWIST where since x = {p, cos(#))}, K is a function of 4 variables 
p, cos(#), p' , cos(#'). It is difficult to estimate a 4-dimensional function accurately 
from Monte-Carlo. In addition, a general feature of unfolding methods is a need for 
a regularization parameter, which biases the result. 

Another approach is to approximate (|7.ip by an analytic expression, and use 
the resulting / rec (x;A) to fit the data. Some terms in the approximate expression 
usually need to be determined from simulations. This method is used by e.g. [27]. It 




(7.1) 
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would be difficult to find a suitable expression for the high precision representation 
of the spectrum required by TWIST. 

The rest of this chapter describes yet another technique, which was used by 
TWIST. The idea is to parametrize the reconstructed spectrum in terms of A, and 
use that parameterization to fit the data. The method is similar to the one described 
in [83]. The TWIST method differs from [83] in the fact that we fit only spectrum 
shape parameters, but not the absolute normalization. 

7.1 The fitting method 

Appendix 110.41 derives an expansion of a reconstructed spectrum in the general 
case. In TWIST background contamination is < 10~ 4 [S3]. Moreover, the primary 
background is simulated by TWIST Monte-Carlo. This means that there is a can- 
cellation of the corresponding spectrum distortion, and the effect is negligible at the 
10 -3 level of precision. Therefore Eq. ()C.22p can be simplified to 



rii(X + AA) 



(A) + X>A Q £-V (7.2) 



7l i> 

a=l 



a=l 

where we also have omitted the 0(AA 2 ) term. Here rij(A + AA) are bin contents of 
a normalized data histogram, nj(A) are corresponding "base" Monte-Carlo values, 
vf are bins of a Monte-Carlo histogram for a "derivative" spectrum corresponding 
to A a , and £~ l is a constant that can be determined from simulation. See Appendix 
[TO for details. 

To extract values of the muon decay parameters, we minimize the 

X 2 = £ ( "' } (7-3) 



n 



?? 



MC 



where n,p ata are the normalized contents of a data spectrum histogram, and 
are calculated according to (17.20 using Monte-Carlo "base" and "derivative" spectra. 
The errors are assumed to be Gaussian, since the available statistics are sufficiently 
high 1 . The error on the content of bin i of an input histogram is taken as the square 
root of the number of entries in that bin, then the errors on different MC histograms 
are combined following (|7.2p . and of is calculated as of = (o-f ata ) 2 + (of 10 ) 2 . 



The statistics are always high for "data" and "base" spectra. Some derivative functions cross 
zero in the fiducial, and around the crossings the count of events in a bin for the corresponding 
derivative may be small. However the total error on the bin is dominated by "data" and "base" 
distributions, so the smallness of a derivative count has no importance. 
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The technique of extracting muon decay parameters through expansion (|7.2p 
takes into account effects of the detector (such as interactions of decay positrons 
in the detector materials) and of the reconstruction (such as possible biases and 
inefficiencies of the track fitting). No further corrections to the fit result are required. 

An important advantage of the method is that effects of the reconstruction can- 
cel exactly, since the same software is used to process real data and Monte-Carlo. 
This fact allows for attributing all systematic effects to deficiencies of the simulation. 
(Indeed, a perfect simulation would be reconstructed exactly as data by any recon- 
struction software.) On the other hand, the TWIST detector is very thin, so that 
the spectrum distortions it causes are small in the first place. Thus any discrepancy 
of the simulation is multiplied by a small factor, and that helps to achieve a high 
precision on the final result. 

The fitting technique assumes independently reconstructed decays (the limit of 
zero beam intensity). If there is more than one muon decay in an event, they are not 
reconstructed independently, and, strictly speaking, the MC integration formulas 
become invalid. Another interpretation is that reconstruction of one decay from an 
event is affected by the presence of the second, so the response function K becomes 
dependent on A. A systematic error accounting for beam intensity effects has been 
evaluated, see chapter 

It is beneficial to choose parameters A so that F is linear in A. Then does 
not depend on A, and the same Monte-Carlo derivative spectra can be used for any 
base A. Here Ni is the number of entries in bin i of the spectrum histogram, see 
Appendix 110.41 This is why instead of A = {p,P^,5} TWIST used A = {p, z,w}, 
where z = P^£,\p^s=const, and w = P^S. However even with the linear parametriza- 
tion = the expansion (|7.2p is still an approximation, because -q^- 7^ if the 
normalization is affected by A. This bias could be overcome by doing iterations: 
generate MC spectra using A = X^ n \ do the fit, take the fit result as A^ n+1 ) and 
repeat. In practice Michel parameters p, £, 5, were all already known to better than 
10~ 2 precisions before the TWIST measurement, so that the 0(X 2 ) contribution 
could be brought down below the 10 -4 level in the first fit, and no more iterations 
were required. 
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It is easy to show that the second order term in the expansion is 

AX a AX/3 f 1 d 2 N t N t d 2 N 
2 \NdX a d\p ~ WdXadXp 

1 dN dNi 1 dN dNi 2N { dN dN 1 
~ N^dKdX^ ~ N^dX^dX^ + WdX^dX^j 

The first two terms in the braces contain a second derivative of the rate and vanish 
in the case of a linear parameterization. The three remaining terms all contain a first 
derivative of the total number of reconstructed events. In the case of a symmetric re- 
sponse function K(E, cos(0), E', cos(#)') = K(E, — cos((9), E' , — cos(#)') derivatives 
dN/d(P^\ P ^ s ) = dN/d(P^5) = 0, so that in fact the dependence ([12]) of n { on 
P[i£,\p^5 and P^S is exact without the 0(AX 2 ) term, and the only deviation comes 
from (Ap) 2 . 

7.2 Specifics of 5 

Two ways of extracting 5 from data were used by TWIST. A 2-dimensional spectrum 
(|1.3p provides the most detailed information, and can be fit to extract p, £, and 5 
simultaneously. (The rj parameter was fixed to the world average, because TWIST 
could not provide a competitive constraint on it.) All systematics were evaluated, 
and the final result extracted, using this approach. 

From (|1.3j) - (|1.5j) it is clear that an "upstream minus downstream" spectrum, 

/(p,cos(0)) - /(p,-cob(0)) oc F AS (p) cos(e) (7.5) 

is manifestly independent of p and rj. Of course, for a reconstructed spectrum it is 
true only to the degree that the response function is symmetric, 

K(p,cos(9),p',cos{6')) = K(p,-cos(9),p',-cos{6')). (7.6) 

Integrating out cos(#) in (17. 5p over a fiducial region, we obtain a 1-dimensional 
momentum spectrum that can be fit with only 2 parameters, A = {P^\p^S: PfiC^}- 
Such fits were employed to cross-check the result, chapter [9) 

The previous measurement of 5 [27] fitted instead the momentum dependence of 
the asymmetry A(p) oc F^{p)/Fis(p). However A{p) is a 1-dimensional spectrum, 
which still depends on all the 4 muon decay parameters. A value of p had to be 
assumed to extract 6 in [27]. TWIST did not use asymmetry fits because of the 
disadvantage of this method. 
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Obviously 5 could not be measured with unpolarized muons, because it only 
enters in the angle-dependent part of Eq. (I1.3p . However, the knowledge of the 
absolute value of muon polarization P^ was not required to determine 5. Any 
inadequacies of simulating muon depolarization in the beam line and the detector 
could be absorbed into P^£, which is a free fit parameter. To achieve this, we need 
to rewrite (jl.5p replacing 



TWIST uses radiative corrections computed within the Standard Model, where 
£ = 1, therefore (|7.7p does not introduce new assumptions. However after this 
replacement Eq. (jl.3p contains only a single parameter P^, instead of separate P^ 
and £. 

7.3 Tests of the fitter 

TWIST implementation of the fitting technique (17.2p ~ (17.3p is based on the ROOT 
rewrite [75] of the MINUIT package [85] . 

To test the program, about 3 x 10 11 muon decays were sampled for "data" spectra 
using p = 0.76, ij = 0, P^ = 1,5 = 0.76, the same amount for "base" spectra using 
p = 0.74, r] = 0, P^ = 0.97, 5 = 0.74, and 10% of that amount for each of the 
p, i], P^\p^s, P^S derivatives. (The derivative spectra do not depend on Michel 
parameters.) Then different fits were done using these decays. The fiducial region 
for the tests was defined as 



All histograms in the tests had 110 bins in momentum from to 55MeV/c, and 100 
bins in cos(#) from -1 to 1. This is the same binning as used in the actual data 
analysis. 

For one test, the decays were split into samples of equal size. The size of the 
sample, 4.8 x 10 7 decays, was chosen as to obtain approximately 10 7 "data" events 
in the fiducial. Each of the "data" samples was fit to a different "base" sample using 
a new set of "derivatives", so that results of all the fits are statistically independent. 
Tests were made in which all 4 parameters were fit as well as tests in which r/ was set 
to zero so that only 3 parameters were fit. The fits were performed in the P^p^s, 
P^5 parametrization and the results were converted to the P^, 5 parametrization 




(7.7) 



20MeV/c < p < 50MeV/c 



0.54 < |cos(0)| < 0.80. 



(7.8) 
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Figure 7.1: Distributions of: (top) fit probability, (middle) deviation of the param- 
eter 5 from the true value, (bottom) reported fit error on 5. Each entry in the 
histograms comes from a statistically independent fit using the same number of 
events. 
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Parameter Mean bias x 10 4 Mean/ (Error on mean) 



P 



8 



-0.026 
-0.014 
+0.287 



-0.24 
-0.06 
+1.45 



Table 7.1: Absolute and relative biases for different fit parameters. 

using formulas from Appendix 110.41 The distribution of fit probability (computed 
from fit x 2 an d the number of degrees of freedom) for the case of the 3 parameter 
fits is shown on the top plot of Fig. 17.11 It is flat, as expected. The biases on all 
fit parameters are consistent with zero, see table 17.11 and the estimates of fit errors 
are close to the RMS of the corresponding distribution. As an example, the middle 
and the bottom plots of Fig. 17.11 show distributions of 5 fit — <5 truc , and of the error 
as, respectively. 

Another test looked at the performance of the fitter as a function of statistics. 
The size of a "data" sample was varied from 10 6 to 2 13 x 10 6 events. (About 2 x 10 5 
to 2 x 10 9 events in the fiducial.) For each sample size, 18 fits were performed 
using the same size of the "base" sample, and 10% of that size for each of the four 
derivatives. Each point on Fig. I7.2H7.4I aggregates 18 fits. Again, all of the fits in 
this test were statistically independent. Fig. 17.21 demonstrates that fit errors, except 
for the lowest tested statistics, scale as l/y/~N. No statistically significant biases 
were observed in the test, see Fig. 17.31 It can be seen from Fig. 17.41 that fit errors 
are underestimated when the statistics is low. However they are consistent with the 
spread of the fitted parameters when the statistics used is higher than about 10 6 
events in the fiducial region. Our measurement used more than 10 data events per 
typical fit, the lowest statistics fit had 0.79 x 10 events in its fiducial region. Thus 
the fitting technique (|T.2|) — (jT.3[) . the conversion formulas from Appendix 110.41 and 
the software implementation of the fitter, have been completely validated for the 
measurement. 

7.4 Blind analysis 

Blind analysis is an increasingly popular tool to avoid (subconscious) experimenter's 
bias when doing a physics measurement. There are subjective decisions to be made 
in e.g. setting the cuts and rejecting "bad" data samples. Several different choices 
may be equally valid and what gets actually used may be affected by the knowledge of 
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Figure 7.2: Scaling of fit errors with statistics. Horizontal axis: the number of data 
events in the fiducial region. Vertical axis: the mean value of reported fit error of 
18 fits at the given statistics. Top to bottom: p, rj, P^£, 5. 
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Figure 7.3: Normalized biases as functions of statistics. Horizontal axis: the number 
of data events in the fiducial region. Vertical axis: (Mean)/(Error on mean) of the 
18 fits at the given statistics. Top to bottom: p, rj, P^£, 5. 
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Figure 7.4: The reliability of fit error estimate as function of statistics. Hor 
izontal axis: the number of data events in the fiducial region. Vertical axis 
RMS/(Mean fit error) of the 18 fits at the given statistics. Top to bottom: p 



53 



1 



0.75 



0.5 
(0 

> 0.25[ 



1950 1955 



1960 

Year 



1965 1970 



Figure 7.5: Experimental determination of the Michel parameter p since 1950. The 
solid line represents the V — A value p = 3/4, 

what version gives a better agreement with the expected answer. Another possible 
source of a bias is looking for software bugs, or additional sources of systematic 
uncertainty when a result does not agree with the expectation, and not looking for 
them as hard otherwise. Often considerable judgment is involved in estimating the 
size of systematic uncertainties. Knowing how close a result is to an expected answer 
may affect the quoted error. A good discussion of motivation for blind analysis, and 
more examples, can be found in [86J. 

There is evidence for such bias in some particle physics measurements. For ex- 
ample, "history plots" in [87] show non-statistical variations of several measured 
quantities with time. In [23], there is the following remark about history of mea- 
surements of the Michel parameter p, which is shown on Fig. 17. 5t 

The curve shows the improvement of the experiments, but perhaps also 
the prejudice of the experimentalists. 

The point is that human bias may introduce an unquantifiable systematic uncer- 
tainty in the result of a measurement. It is desirable to avoid the possibility of 
such a bias. This can be accomplished by doing analysis in a "blind" fashion, i.e. 
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Figure 7.6: TWIST blind analysis scheme. 

keeping the final result hidden till the analysis is essentially complete. The value 
of a measurement does not contain any information about its correctness and is of 
no use in performing the analysis, therefore hiding the answer does not impede the 
work. 

TWIST implementation of a blind analysis 

Among our requirements for a blind analysis scheme were: 

• Does not exclude any TWIST member from doing any part of the analysis. 

• Convenient to use. 

• Minimal modifications to the existing software. 

• Hard to break. 

A scheme of implementation satisfying these criteria is shown on Fig. 17.61 The idea 
of the method is to blind the MC samples, not the fitter. It is clear from (|7.2p that 
the fitting method gives only deviations of the Michel parameters in data from the 
values used to generate a base Monte-Carlo spectrum. Thus it was sufficient to hide 
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the values of Michel parameters used for MC production. The secrecy was based 
on using an asymmetric cryptography. A private-public key pair was produced, and 
the private key locked up in a place not accessible by TWIST members. 

The piece of software written specifically to make a blind analysis possible is 
micheld, which is essentially a muon decay spectrum generator. The spectrum it 
produces includes all radiative corrections described in chapter |U The program runs 
on a computer which is not controlled by the TWIST group, and none of TWIST 
members could login there during the analysis period. (Symbolized by a "wall" 
around micheld on Fig. 17.61 ) micheld is a multi-threaded TCP/IP server, accepting 
data and control requests from a network. An operator uses micheld_ctl to instruct 
micheld to produce a random set of Michel parameters. They are sampled uniformly 
within the following limits: 

p=^±0.02, 77 = 0, P M £ = 1±0.03, 6 = ^±0.02 (7.9) 

A candidate set of values is tested for being physically allowed (the end point asym- 
metry Pn^S/ p < 1) before it is accepted. An accepted set of parameters is encrypted 
using the public key, and the encryption result is stored in a database. By another 
operator request an accepted set of Michel parameters is used to generate a series 
of muon decay samples, which are written to disk. 

During a Monte-Carlo production run, a GEANT process obtains a sample of muon 
decays from the disk through micheld. Every time GEANT needs to decay a muon, it 
uses the energy and angle (with respect to the spin of the muon) of the next decay 
in the sample. Since different muon decay samples were produced with the same 
(unknown) Michel parameter values, we had the possibility to study consistency 
between different data sets, and to estimate systematics by fitting one MC sample 
to another, as explained in chapter [SJ 

After the analysis was complete, the "black box" was opened, and the final values 
of parameters were computed using the results of the fits and the revealed MC val- 
ues. A small complication arises from the fact that the (-P^£|p £5, P^S) — > (P M £,<!>) 
conversion (Appendix ll0.4p requires the knowledge of "true" MC parameters. That 
was addressed by using the known approximate numbers P^o = 1, 8q = 0.75, at 
the "blind" stage. This approximation introduces an uncertainty of about 5% on 
the deviation A5 from the fit, which translates into a 5% uncertainty on sensitivi- 
ties of 5 to different systematics (see chapter [8]) . This was adequate for doing the 
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analysis. After opening the black box all data and systematics fits were re-run using 
the revealed values of P^o an d Sq to get rid of the additional uncertainty. This 
was a mechanical procedure not involving any judgment, so it did not violate the 
philosophy of blind analysis. 
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Chapter 8 

Determination of systematic 
uncertainties 

To estimate the systematic uncertainty of the result, it is essential to account for 
all sources of systematic effects. On the other hand, it is important to avoid double 
counting, so that the same physical cause of a bias is not included more than once 
in the estimate. Also, one has to distinguish evaluation of systematic uncertainties 
from consistency checks [88]. The decision on what effects to consider and exactly 
how to treat them involves judgment and is, to some degree, arbitrary. An important 
feature of the present measurement is that it was done using a blind analysis scheme. 
(See l7.4i ) This means that a complete list of systematic effects, along with a method 
to evaluate each of them, was fixed before "opening the black box" and revealing 
the measured value of 5. Such an approach reduces the possibility that the obtained 
estimate of systematic uncertainty is subjectively biased. 

In our approach (chapter [7]) all systematic effects can be attributed to imper- 
fections of the Monte-Carlo. A simulation perfectly reproducing data would be 
reconstructed exactly as data by any reconstruction software, thus the result of a fit 
of data to Monte-Carlo would be unbiased. In other words, effects of reconstruction 
cancel in the comparison of data to Monte-Carlo to the degree that the simulation 
reproduces the data. Therefore the list of systematics does not include effects of the 
reconstruction. This of course does not mean that the quality of the reconstruction 
software is irrelevant: the sensitivity of the result to a given imperfection of the 
simulation may be reduced by improving the reconstruction program. 

For TWIST the possible sources of systematic uncertainties can be classified into 
the following independent groups: positron interactions, spectrometer alignment, 
chamber response, momentum calibration, and muon beam stability. Some of the 
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specific effects within these groups bias the result in the same way for all the data 
sets, while other effects are time dependent and may contribute differently to the 
different data sets. The former effects only contribute to the uncertainty of the 
result, while the time-dependent uncertainties also should be used in computing the 
weighted average value of 5 from the four data sets, (chapter 0) 

Most of the systematics effects were evaluated using the following technique: 
evaluate the sensitivity R' = d5/da of the result 5 to a systematic parameter a known 
to a precision ±o" a , and add R'a a quadratically to the systematic uncertainty [88] , 
To estimate R', a data set set was taken or a Monte-Carlo set was generated under 
a different condition a tcs t 7^ ^nominal- The reconstructed "test" spectrum is then fit 
to the "nominal" one using (j7.2f ) (|7.3f) . This expresses the change in the spectrum 
shape due to the systematic effect in terms of changes in the Michel parameters. 
The systematics estimate can therefore be written as 

R'ffa = ~T — &a = (<5test ~ ^nominal) 7 T = 7; (^test — ^nominal) (8-1) 

£± a (.Qtest — a nominalJ >-> 

where we have introduced the scaling factor S = (atest — «nominai)/fa- 

The following subsections present in tabular format summaries of individual 
systematics for each group, followed by a short explanation of each entry in the 
table. If all systematic uncertainties R'a a are identical for all data sets used in the 
measurement, a single column is used to present them, as in Tables 18.11 18.21 18.51 
Otherwise individual numbers for data sets A, B, 1.96 T, and 2.04 T are shown, 
Tables E3J El 

8.1 Spectrometer alignment 



Table 18.11 summarizes alignment related systematic uncertainties, all of which are 
data set independent because the TWIST detector was mechanically stable and its 
parts did not move during the data taking [52] . 



Name 


10 ;i x AS 


Scaling 


10 a x R'a a 


Translations 


0.39 


28 


0.01 


Rotations 


-4.33 


39 


-0.11 


Z (longitudinal) 


-1.07 


10 


-0.11 


B field to detector axis 


-1.86 


3.1 


-0.60 


Total 






0.62 



Table 8.1: Alignment systematics. 
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The translational alignment of the wire planes relative to each other is measured 
using straight tracks produced by 120 MeV pions with the solenoid off. The 
accuracy of the resulting alignment, a a , is 5 /xm [52]. The test spectrum 
was produced by analyzing the nominal data set B with a special alignment 
file, which was produced by applying random shifts to the nominal alignment 
corrections file. The resulting translational spread of the wire planes from 
their nominal positions was 140 /jm (RMS), giving the scaling factor S = 
140 /xm/5 fim = 28. 

The rotational alignment systematic was determined in a similar way, using a 
specially prepared rotational corrections file. The introduced angular spread 
of 0.39° (RMS) yields a scaling factor of 39 compared to the 0.01° precision of 
the nominal rotational alignment |52j . 

The Z (longitudinal) alignment of the wire planes is estimated to be accurate to 
30^m from mechanical precision of the detector construction [52]. A Monte- 
Carlo set was generated with Z positions of the planes offset by 300 fim (RMS) 
and compared to a nominal MC set, thus producing a scaling factor of 10. 

B field to detector axis. The nominal Monte-Carlo generation and data analysis 
assume a perfect alignment of the detector axis to the coordinate system of 
the magnetic field map. To produce the test spectrum, the field map was 
rotated in GEANT by 0.25°. The actual misalignment is estimated from data 
by fitting (an approximation of) a helix that is not aligned to the detector 
axis to the positron tracks, with the alignment angles being two additional 
free parameters in the fit. The average misalignment found in this way is 
0.08°, so the scaling factor is 0.25°/0.08° = 3.1. 

8.2 Positron interactions 

Energy smearing. This systematic accounts for any mismatch between data and 
MC in the momentum resolution. This mismatch has been observed at the 
end point (chapter [6]) . A test spectrum was produced by applying a Gaussian 
smearing (cr ptSme ar = 200 keV/c) to the transverse component of momentum, 
Pt, of reconstructed Monte-Carlo events. The same Monte-Carlo data analyzed 
in the standard way, without any smearing, gave the reference spectrum. 
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Name 


10 3 x AS 


Scaling 


10 3 x R'a a 


Energy smearing 


0.58 


4 


0.15 


Multiple scattering 


0.10 


20 


0.00 


Hard and intermediate interactions 






0.53 


Detector materials 


-0.73 


2 


-0.36 


Outside materials 


-2.09 


70 


-0.03 


Total 






0.66 



Table 8.2: Positron interactions systematics. 



For o"p iSmear = 50 keV/c the width of the end point, as determined by the 
energy calibration procedure (chapter [6|) , agreed with data, o"g I Q eared MC ~ 
(jgg ta . Therefore the scaling factor used was 200/50=4. 

Multiple scattering. The multiple scattering systematics addresses a potential 
deficiency in the simulation of multiple scattering (chapter 2]). The angle 9 of 
reconstructed Monte-Carlo events was smeared with a Gaussian using 



p(MeV/c) 



3.2 



|cos(0)| 

to produce a test spectrum. Here K is a parameter, while the functional 
dependence on momentum and angle comes from a simplified formula for the 
multiple scattering angle of a relativistic particle in matter. (See e.g. [87], page 



245.) The 



cos(0) 



term is proportional to the amount of material traversed by 



a particle in the planar TWIST geometry. For K = 1 rad, p = 30 MeV/c, 
and cos(0) = 0.7, Eq. JS2J) gives d0 smcar ~ 29 mrad. The size of the discrep- 
ancy between data and Monte-Carlo was estimated as the bigger between the 
differences in mean and RMS between data and Monte-Carlo in the valida- 
tion studies (chapter [4]) . None of the differences exceeded 1.5 mrad [89], so a 
scaling factor of 20 was chosen 1 . 

Hard and intermediate interactions. The systematic uncertainty due to the 
imperfect simulation of hard and intermediate interactions was estimated in 
the following way [90] . A spectrum of reconstructed and thrown positron mo- 
menta was prepared for events that lost less than 1 MeV/c in the detector 
(according to Monte-Carlo information). A similar spectrum was generated 



x This estimate of the scaling factor is not well justified. However the typical smearing angle of 
29 mrad is larger then the angular resolution (chapter [SJ, and the "raw" effect of the smearing is 
still small. This is why the obtained estimate of this systematics was never refined. 
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^omentum^reconsm^ 
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Figure 8.1: Ratios of (reconstructed) /(thrown) momentum spectra for: Apyic < 
1 MeV/c (left), Ap MC > 1 MeV/c (right). Only events with 0.7 < | cos(0)| < 0.84 
are included in the spectra. (Plots from |90j.) 

for positrons that lost greater than 1 MeV/c. The 1 MeV/c number is an ar- 
bitrary boundary between "intermediate" and "hard" interactions, consistent 
with the boundary used in the Monte-Carlo validations (chapter [5]) . Ratios 
of the (reconstructed)/(thrown) distributions, presented on Fig. 18.1} show the 
distortions of the momentum spectra for the two classes of events. 

The fractional yield changes over the range 25-50 MeV/c are 

Us = (si — S2)/norm « 0.0067 (left plot, intermediate interactions), (8.3) 
Vh = (hi — /i2)/norm « 0.0070 (right plot, hard interactions). (8-4) 

Here hi » 0.7458, h 2 ~ 0.7405, si » 0.0153, s 2 ~ 0.0102 are the readings at 
25 MeV/c and 50 MeV/c from the plots, and norm = ^(/ii + /i2 + si + S2)- The 
change of the Michel parameter p by 0.0010 leads to a fractional yield change 
of 0.0018 over the same range. GEANT has been validated to 14% for hard 
interactions and to 5% for intermediate interactions (chapter 2]) . Therefore 
an uncertainty on p can be estimated as jyggjf x (0.05 x y s + 0.14 x y^j- 
To determine the effect on 6, we scale the uncertainty on p by the ratio of 
AS/Ap ~ 0.719 obtained in the "Detector materials" systematic below. The 
final number is 

0.0010 



0.719 x 



0.0018 



x (0.05 x y a + 0.14 x y h ) w 0.00053 



(8.5) 



Detector materials. The nominal thickness of the graphite coating on each side 
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of the stopping target is 10 /im. A test Monte-Carlo set was generated with 
30 /iin graphite coating and fit to a nominal 10 /im MC set. The setting of 
the gas absorber in the test MC run was tuned to keep the muon stopping 
distribution centered in the target, because the systematic effect related to 
shifts of this distribution is accounted for separately. Since the thickness of 
the graphite coating on each side of the stopping target is known to be between 
5 /im and 20 /im [52], a scaling factor of (30 — 10)/(20 — 10) = 2 has been 
applied. 

Another output from this systematic study is the relative size of effects on 
5 and p due to interactions of the positrons in detector materials, AS/Ap ~ 
0.000726/0.001010 ~ 0.719, used above in the estimation of the uncertainty 
due to hard and intermediate interactions. 

Outside materials. Decay positrons after leaving the tracking volume may scatter 
off the outside structures of the detector, re-enter the tracking volume, and 
produce more hits, and consequently confuse track reconstruction. The biggest 
source of backscatters is the upstream beam package, which holds the trigger 
scintillator and the degraders (Fig. I2.ip . There were no materials other than 
air at the downstream end during the normal data taking. 

To estimate the effect of an imperfect Monte-Carlo simulation of the positron 
backscattering process, a special data set was taken with an aluminum plate 
mounted outside of the downstream end of the detector. A fit of this data 
set to a nominal data set produced the shift Ad = —2.09 x 10 -3 , shown in 
Table E2J 

Backscatters and beam particles overlapping in time with a decay positron 
may not be distinguishable on an event by event basis. A study of backscat- 
ter rates in data and Monte-Carlo [91] used the "PC time of flight" variable, 
Tpc = t u — td, where t u and td are the average times of hits in the 4 most 
upstream and most downstream PCs. Accidental overlaps, such as those with 
beam positrons, produce a flat background, while backscatters produce a peak 
in the Tpc distribution. This difference provides a way to measure the rate 
of backscatters. The rate of backscatters from the downstream direction un- 
der nominal conditions was demonstrated to be ~ 0, and the scaling factor 
was estimated as the ratio of the backscatter rate from the downstream alu- 
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minum to the difference of the backscatter rates between the nominal data 
and simulation [91] : 

S = ,„ NDdD Z° ■ - 70. (8.6) 



\N GuS - AW | 

Here NquS is the number of backscatters seen in GEANT from the upstream 
material under the "standard" conditions, Ne> u s is the number of backscatters 
seen in data from the upstream material under the "standard" conditions, 
and NjjdD is the number of backscatters seen in data from the downstream 
aluminum. 



8.3 Chamber response 



Name 


10 s x A5 


Scaling 


10 3 x R'a a 


A 


B 


1.96 


2.04 


DC efficiency 


0.27 


50 


0.01 


0.01 


0.01 


0.01 


PC efficiency 


0.07 


50 


0.00 


0.00 


0.00 


0.00 


Dead zone PC 


0.46 


6 


0.08 


0.08 


0.08 


0.08 


Dead zone DC 


1.38 


15 


0.09 


0.09 


0.09 


0.09 


Up-down differences 


-0.19 


4 


-0.05 


-0.05 


-0.05 


-0.05 


HV variations 


0.08 


20 


0.00 


0.00 


0.00 


0.00 


Temp, Pressure 


-2.66 




-0.35 


-0.35 


-0.22 


-0.35 


Foil bulges 


-1.3 




-0.52 


-0.26 


-0.52 


-0.26 


Crosstalk 


0.01 


10 


0.00 


0.00 


0.00 


0.00 


TO variations 


-1.83 


10 


-0.18 


-0.18 


-0.18 


-0.18 


Total 


0.67 


0.49 


0.61 


0.49 



Table 8.3: Chamber response systematics. 



DC efficiency. A special analysis of the nominal data set B deleted 5% of DC 
hits before passing events through the standard reconstruction chain. Once 
a hit was marked for deletion, all hits on the same wire within 700 ns were 
also deleted, since they were likely to come from the same track. The obtained 
Michel spectrum was fit against the standard analysis of the same set. The 5% 
inefficiency that was introduced for this test corresponds to an exaggeration 
factor of 50, because the actual efficiency of the DCs is about 99.9% [52] , 

PC efficiency. The effect of inefficiencies in the PCs was estimated in a similar way. 
The test spectrum was produced from set B with a 5% artificial inefficiency 
and a hit removal time interval of 50 ns. This systematics becomes negligible 
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after applying the scaling factor of 50. Note that even the "raw" effect of 5% 
PC inefficiency is small, because PCs are not used in the track fitting. 

Dead zone PC. Muons slow down in the detector and become highly ionizing 
before they stop. The large space charge they create in wire chambers may 
"deaden" a section of wire for some time after a passage of a muon. The effect 
is largest in PC 6, the chamber which is closest to the muon stopping target on 
the upstream side, and is also observable in PC 5. The nominal Monte-Carlo 
did not simulate this effect. 

In a special study [92j [93] , the dead zone effect was simulated by introducing 
a 100% inefficiency along the wire around the point of the muon hit. The 
rectangular in space inefficiency zone exponentially shrank in time. The initial 
length of the zone was computed as the length of the projection of the muon 
created space charge on the wire plus AL. A "realistic" simulation with the 
parameters AL/2 = 0.24 cm and Th ca i = 2444 ns reproduced the "dip" in PC 6 
efficiency around a muon hit that was observed in data reasonably well. An 
"exaggerated" simulation used AL/2 = 5.00 cm and Th ca i = 3500 ns. Both 
simulations used the same dead zone parameters for all PC planes. 

To measure the sensitivity R', the "exaggerated" MC was fit against a nominal 
(no dead zone) simulation. The scaling factor was estimated using the ratio of 
the number of PC positron hits lost due to the dead zone in the "exaggerated" 
MC to the number lost in the "realistic" MC. 

Dead zone DC. A similar "dead zone" effect is also expected in the DCs. Its 
magnitude is smaller then in PC 6, because muons are less ionizing farther 
from their stopping position. A DC dead zone is also seen at a smaller solid 
angle from the stopping target than PC 5 or 6 dead zone, making the effect 
harder to observe. 

With no estimate of DC dead zone parameters available from data, the same 
AL/2 and 7h ea i parameters as used in the PC study were applied to the DCs. 
The scaling factor was estimated using the ratio of the number of DC positron 
hits lost due to the dead zone in the "exaggerated" MC to the number lost in 
the "realistic" MC. 

Upstream-downstream differences. The average number of degrees of freedom 
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of the track fits, (ndof), related to the number of drift chamber hits used 
in the tracking, is different between data and Monte-Carlo. Among many 
possible explanations (like the cathode foil bulges, see below) is an inaccuracy 
of simulation of the drift chamber response in the corners of a drift cell. Special 
analyses were run, which excluded "cell corner hits" with t > t max from the 
final track fitting. For a special analysis of Monte-Carlo, i max = 400 ns was 
used. In the special data analysis, the cut was tuned to match the average 
(ndof) of the data track fitting to that of the 400 ns Monte-Carlo analysis. 
The tuning resulted in t max = 522 ns. 

The data spectrum produced in the special analysis was fit to a nominal analy- 
sis of the same set, giving AS = 0.000 x 10 -3 . A similar fit of special to nominal 
Monte-Carlo spectra gave AS = 0.193 x 10~ 3 . The "raw" effect was taken as 
the difference between the data and the Monte-Carlo fits: AS = —0.193 x 10~ 3 . 

The effect of corner cell inefficiencies on the spectrum could be taken with a 
scaling factor of 1. However while tuning the average (ndof) = \ ((ndof) + 
(ndof) dre ), the long drift time cuts exaggerated a difference in the asymmetry 
^ndof = ((ndof) up - (ndof) dn )/( (ndof) up + (ndof) dn ) between data and Monte- 
Carlo by about a factor of 15, from A^f - ~ -0.16% for the nominal 
analyses to 2-3% for the long drift time cut analyses. Since S is an asymmetry 
parameter, a scaling factor of 15 was another possible choice. It has been 
decided to use a factor of 4 (the geometric mean) for this systematic. 

HV variations. This systematic represents the effect of our imperfect knowledge 
of high voltage on the wire chambers. The test spectrum was produced by 
analyzing the nominal data set B using drift tables corresponding to 1850 V, 
and comparing it to the same data set analyzed with the nominal, 1950 V, 
drift tables. The scaling factor is 20 because the accuracy of the high voltage 
if 5 V. 

Temperature and pressure. This systematic uncertainty represents effects of vary- 
ing gas density in the TWIST drift chambers, caused by variations of the 
atmospheric pressure and outside temperature. 

A special Monte-Carlo set was generated with settings corresponding to the 
temperature of —10° C, instead of the nominal +20° C, and fit to a nominal 
Monte-Carlo set, yielding AS = —2.66 x 10~ 3 . Scale factors were determined 
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individually for each data set by comparing variations in the gas density that 
occurred during data taking (available from the Slow Control DAQ records) 
to the « 10% density change in the —10° C test MC set. The scaling factors 
are: 10 (set A), 10 (set B), 20 (1.96 T), 10 (2.04 T). 

A change in the gas density also affects the muon stopping distribution. To 
avoid double counting the stopping distribution effect, the latter systematic 
number, 0.09 (discussed below), was linearly subtracted from scaled per-set 
temperature and pressure estimates. 

Foil bulges. During the data taking, the differential pressure between the drift 
chambers and the enclosing He/N2 volume was not always stable. That led 
to a movement of the cathode foils, affecting both the space-time relationship 
of the drift chambers, and the average number of wires in a plane hit by a 
track at a given angle. For example, in the nominal geometry with square 
4 mm x 4 mm drift cells a track at 6 < 45° can never hit more than 2 cells in 
a plane. But if the cathode foils bulge out, extending the cell size along the 
detector axis, the same track might produce more hits per plane. 

To estimate the systematic uncertainty due to the foil movements, a special 
Monte-Carlo set with the foils moved out by 500 fim was produced and fit to 
a nominal MC set, giving A5 = —0.00130. 

Several measurement of the differential pressure, done during the running 
period, were correlated to an analysis variable derived from the most probable 
value of x 2 / n dof in the tracking and corrected for gas density effects |94j . That 
variable was then used to set the following scaling factors for different data 
sets: 2.5 (set A), 5 (set B), 2.5 (1.96 T), 5 (2.04 T). 

Crosstalk. The test spectrum was produced by turning off the crosstalk removal 
algorithm (page [26]) and re-analyzing nominal set B. As expected, the effect is 
very small since decay positrons are weakly ionizing and do not produce large 
signals in the chambers, which could induce crosstalk. The removal algorithm 
is estimated to be correct at least 90% of time, thus the scaling factor of 10. 

TO variations. Time offsets TO for different wires were calibrated using 120 MeV 
pion tracks with solenoid off. Results from a calibration run taken at the 
beginning of the data taking period were used in the nominal data analysis. 
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To estimate the systematic uncertainty associated with the calibrations, an- 
other TO run was taken at the end of the data taking period. The per-channel 
differences in the extracted offsets between the two runs were multiplied by 
10 (the scaling factor) and added to the nominal TO values. The obtained TO 
calibration file was used to analyze data set B to produce the test spectrum, 
which was fit against the nominal set B spectrum, resulting in A5 = —0.00183. 

8.4 Momentum calibration 



Name 


10 3 x A5 


10 3 x R'a a 


A 


B 


1.96 


2.04 


End point fits 




0.27 


0.21 


0.21 


0.24 


Field map 


0.68 


0.07 


0.07 


0.17 


0.34 


Total 


0.28 


0.22 


0.27 


0.42 



Table 8.4: Momentum calibration systematics. 



End point fits. Sensitivities of 5 to the energy calibration parameters /3, a u , and 
(chapter [6]) were determined by fitting a nominal spectrum to a spectrum 
with the appropriate energy calibration constant offset by 100 keV/c. Co- 
variance sub-matrices V for (/3,a u ,ctd) from energy calibration fits to data 
and to the corresponding MC spectrum were added, and the resulting matrix 
multiplied with the vector of sensitivities A in the usual way [95] : 

a 2 = A T (V D ^ + V MC )A (8.7) 

Field map. A mismatch between the measured B z component of the nominal 2 T 
spectrometer field and the OPERA field map was fit using the expression 
AB Z = C2Z 2 + C32; 3 + c r r [96]. This simple function describes the residuals to 
within 1 G throughout the entire tracking region, and to within 0.5 G over 
most of the tracking region. Nominal data set B was re-analyzed using a test 
field map, which was prepared by adding 10 x AB Z to the nominal OPERA 
map, and then fit against the standard analysis of the same data. Thus the 
scaling factors are 10 for the 2 T sets A and B. Nominal analyses of the 1.96 T 
and 2.04 T data sets were done using a scaled 2 T field map. Comparisons 
of the scaled versions to the actual B z measurements at those fields gave the 
scaling factors of 2 for the 2.04 T data set, and -4 for the 1.96 T data set. 
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8.5 Muon beam stability 



Systematics related to the stability of muon beam parameters are "data set depen- 
dent" by their nature. Since none contributed significantly to the final result, a 
common estimate for each of the scaling factors based on the worst case data set 
was used, instead of assigning individual scalings to different data sets. 



Name 


10 3 x AS 


Scaling 


10 3 x R'a a 


Stopping location 


0.52 


6 


0.09 


Beam intensity 


0.26 


6 


0.04 


Channel magnets 


-1.29 


50 


-0.03 


Total 






0.10 



Table 8.5: Muon beam systematics. 



Stopping location. The average position of muon stops in the stopping target 
affects the amount of the target material seen by decay positrons. E.g. if 
muons stop before reaching the center of the target, positrons going upstream 
will be less affected by energy loss and multiple scattering than those going 
downstream. The energy calibration procedure (chapter [6]) compensates, to 
first order, for differences in energy loss. The "stopping location" systematics 
covers any remaining effects. 

A special data set was taken with the muon stopping position displaced slightly 
upstream (by introducing more CO2 in the gas degrader). That set was fit 
to a nominal data set to measure the effect. The scaling factor was obtained 
by comparing the ratio adis/a sum for the special set, —0.12, to the spread of 
that ratio for other data sets, ~ 0.02. Here a^s and a sum are the energy 
calibration variables, see chapter [6j 

Beam intensity. The signal rate on the TWIST trigger counter was recorded dur- 
ing the data taking. After rejection of bad runs, the spread of the average 
trigger rate for different runs within the four nominal data sets was found to 
be smaller than 0.6 x 10 3 s _1 . To measure sensitivities of the Michel parame- 
ters to beam intensity, a low rate (1.1 x 10 3 s _1 ) and a high rate (4.7 x 10 3 s^ 1 ) 
data set were taken and fit against each other. The exaggeration factor S in 
this measurement is (4.7 — l.l)/0.6 = 6. 

Channel magnets. This systematic accounts primarily for instabilities in the B2 
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Effect 



Uncertainty 



Spectrometer alignment ±0.00062 

Chamber response(ave) ±0.00056 

Positron interactions ±0.00055 

Stopping target thickness ±0.00037 

Momentum calibration(ave) ±0.00030 

Muon beam stability(ave) ±0.00010 

Theoretical radiative corrections ±0.00010 
Upstream/Downstream differences(ave) ±0.00005 



Total 



±0.00112 



Table 8.6: Contributions to the systematic uncertainty for 5. For set-dependent 
systematics the average values, denoted by (ave), are shown. 



beamline dipole (Fig. l2.2| ). which directly affects the position of the muon beam 
as it enters the TWIST spectrometer. The other dipole magnet, Bl, defines 
the momentum of the muons accepted by the channel, and its instabilities are 
included in the "Stopping location" systematic. 

The strength of the magnetic field in B2 was continuously monitored and was 
stable to 0.2 G. A test data set was taken with B2 intentionally offset from 
the nominal value by 10 G, giving the exaggeration factor of 50. 

The deflection of beam particles by M13 quadrupole magnets is small com- 
pared to the 60° bend by B2. Since the systematic effect of B2 is already 
small, contributions of the quadrupoles to the systematic were neglected. 

8.6 Summary of systematics 

Table 18.61 shows a summary of the systematic uncertainties. The "stopping target 
thickness" is shown separately from the rest of "positron interactions" uncertain- 
ties, and "Upstream/Downstream differences" separately from the rest of "Chamber 
response", following [97] . 

An entry not discussed above is the uncertainty from theoretical radiative cor- 
rections. Theoretical uncertainty on 5 is estimated as 1 x 10~ 4 [18], if terms of up 
to 0(a 2 L 2 ) are included in the spectrum. Here L = l^m^/m 2 ) ~ 10.66, and a is 
the fine structure constant. TWIST uses an even more precise spectrum description 
(chapter [4|) , therefore this estimate provides a safe upper bound on the uncertainty. 



70 



Chapter 9 

Determination of decay 
parameters 

According to the philosophy of a blind analysis, the decision on which data sets to 
use for the extraction of a final result had been made prior to "opening the box" and 
revealing the physics result. Within each of the chosen data sets, an identification 
and exclusion of bad runs was also done at the blind stage of analysis. 

All fits used for the extraction of 5 were done using 2-dimensional histograms 
of reconstructed data and Monte-Carlo spectra in momentum and cos(#), which 
is the complete information available from the detector. The fits were done in the 
linear parametrization {Ap, Az, Aw}, where z = P^S,\p^s=const, an d w = (Sec- 
tion 17. ip . then the results were converted to the usual {p, P^,5} parametrization 
and the covariance matrices recomputed. The 3-parameter fits, unlike the fits to the 
spectrum asymmetry A{p) = -Fas (p) /-^is (p) (see Eqs. I1.3H1.5I) used in the previous 
measurement [27], do not require making any assumption regarding the value of p 
in order to find 5. Sensitivity of <5 to the value of r/ = —0.007 ± 0.013 |87] assumed 
in MC production was checked and found negligible. 

The measurement of the decay parameter 5 uses the following four data sets: 
set A, set B, 1.96 T, and 2.04 T. (See Table ED) Table O shows results of fits 
to the chosen data sets, computed using black box offset values po = 0.74766, 
P M £ = 1.0148, 5 = 0.73645. 

Correlation coefficients for set B are shown in Table. 19.21 Correlations for other 
surface muon sets are very similar. The small (less than 10%) correlation between 
p and 5 confirms that in our approach the two parameters are independent. 

The simulation describes the data well, as it is demonstrated by x 2 /NDOF « 
1 and the reasonable fit probabilities shown in Table 19.11 Figures I9.1H9.5I show 
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Data Set 


5 


P 


x 2 


Probability 


Set A 


0.75087 ± 0.00156 ± 0.00073 


0.75083 ± 0.00083 


1924 


0.27 


Set B 


0.74979 ± 0.00124 ± 0.00055 


0.74911 ± 0.00066 


1880 


0.54 


1.96 T 


0.74918 ± 0.00124 ± 0.00067 


0.74956 ± 0.00066 


1987 


0.05 


2.04 T 


0.74908 ± 0.00132 ± 0.00065 


0.75203 ± 0.00071 


1947 


0.16 



Table 9.1: Fit results. Set-dependent systematic uncertainties from Tables [873]. I8.4| 
18.51 are shown for 5 after the statistical errors. Results for p are consistent with our 
previous measurement |98j . Large depolarizing effects in the graphite coated Mylar 
target (chapter 5|) made the present data unsuitable for an improved measurement 
of P^£, therefore this parameter is not shown. Each fit has 1887 degrees of freedom 
(NDOF). The last column is the fit probability computed from \ 2 an d NDOF. 









AP^ 


A<5 


Ap 


0.157 


0.262 


Ap 


0.157 


0.097 






0.422 


AP^ 




-0.541 



Table 9.2: Correlation coefficients. Left: in the original fit parametrization. Right: 
converted to the usual muon decay parameters, according to formulas from Ap- 
pendix 110.41 

residuals of the fits, providing more details on fit quality. On each of the figures, 
the top left panel shows the normalized deviation of the best fit, (Data — Fit)/<7, for 
each bin of the 2-dimensional fit histogram. The deviation of a bin is normalized to 
its statistical error a, and shown on a color scale. The two solid contours delimit 
the fiducial region (page 1331) . It can be seen that most bins within the fiducial agree 
to 1-2 cr. 

The 2-dimensional fit and data histograms were independently projected onto 
the momentum (top right), and cos(#) (bottom left) axes, ignoring bins outside 
of the fiducial region. Deviations between the obtained data and fit projections 
are shown using the solid marker. The empty marker points on the top right plot 
were obtained by removing the p < 50 MeV/c fiducial cut, and projecting the 2- 
dimensional distributions from the extended region. Similarly the 0.50 < | cos(0)| < 
0.84 cut was removed to obtain points outside of the fiducial on the bottom left 
panel. The bottom right panel shows a match between the data and the fit for the 
Fas{p) distribution, which is the most relevant for 5. It was obtained by angle- 
integrating the spectra separately in the upstream (cos(#) < 0) and the downstream 
(cos(0) > 0) parts of the fiducial region, and computing the difference Fas(p) oc 
Upstream — Downstream. Most of the residuals for bins inside the fiducial region 
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are within 2 a off zero for all the views, and no pattern indicates a systematic 
difference between the data and simulation. 

A series of consistency checks were done. A fit to the Cloud set gave 5 = 
0.75245 ± 0.00526 (stat). The fact that 5 extracted from a data set with an opposite 
(and small) muon polarization, _p^ loud ~ +0.25, is consistent with the value extracted 
from surface muon data, with p™ rface ~ — 1 ; demonstrates the absence of detector 
asymmetries that would lead to different biases on 5 for the two cases. 

Another check involved generating a Monte-Carlo set with values of the muon 
decay parameters determined from set B, and using the produced spectrum to do 
another fit to set B. The fit yielded all deviations of the decay parameters consistent 
with zero, as expected 1 . 

Yet another check used a modified fitting procedure. Instead of doing a 3- 
parameter fit to the 2-dimensional spectrum, a 1-dimensional distribution propor- 
tional to Fas(p) was extracted from it. (See (|1.3p - (|1.5p .) The shape of that distri- 
bution is manifestly independent of p and i] (provided the detector response function 
is symmetric). This 1-dimensional distribution was fit using the P^Ip^s and P^S 
parameters, and an alternative value of 5 was extracted from that fit. A comparison 
of the alternative values, 5 U( j, with the values extracted from fits to the 2-dimensional 
(p, cos(#)) distribution, 52D-, is shown in Table I9T31 The expected variance of the 
difference a^ iS for correlated data was computed as a\ m = |er^ d — ^Id I ' assuming 
that one of the estimators saturates the Minimum Variance Bound [88[ 199] . It can 
be seen from Table [931 that results of the alternative technique are highly consistent 
with those given by the 3-parameter fits. 

We compute the central value of 5 as a weighted average [87] . using for the 
weights a quadratic sum of statistical and set-dependent systematic uncertainties 
from Table 19.11 In the calculation of the final systematic uncertainty we do not 
assume that it shrinks in the combined measurement, and quadratically add set- 
independent and average values of set-dependent systematics as shown in Table 18.61 

lr That test failed for 8 in the first round of TWIST analysis. The problem was traced to a failure 
of implementing Eq. (|7.7[) in the spectrum generator. Since p is decoupled from the asymmetry 
parameters, this flaw had no impact on p, and its value was published [98]. On the other hand, it 
introduced an additional systematic uncertainty on 5, dependent on the difference of the average 
muon depolarization in data and Monte-Carlo. This systematics was estimated to be < 0.001. 
This large value necessitated a reanalysis, with the generator fixed. A new black box was created, 
and a second round of analysis performed, which is presented in this work. Because the effect 
of the mistake on 5 was large, we did not know the value of 5, so the analysis was still blind. 
Other changes between [55] and this analysis are improvements in track selection (chapter [SJ) and 
rotational alignment of the drift chambers. 
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Data Set 


<5ud — <^2D 


Cdiff 


(<5ud — <52D)/Cdifr 


Set A 


0.000218 


0.000293 


0.74 


Set B 


-0.000006 


0.000230 


-0.03 


1.96 T 


0.000168 


0.000228 


0.74 


2.04 T 


-0.000007 


0.000239 


-0.03 


Cloud 


0.000131 


0.000711 


0.18 



Table 9.3: Per-set differences between 5 u( j extracted from a fit to the "upstream 
minus downstream" distribution Fas{p)i an d &2D from a 3-parameter fit to the 
(p, cos(#)) spectrum. 



The final result is 

5 = 0.74964 ± 0.00066 (stat.) ± 0.00112 (syst.) (9.1) 
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Figure 9.1: Fit residuals for set A. Top left: color-coded residuals in the cos(#) 
vs momentum plane. Top right: projection on the momentum axis. Bottom left: 
projection on the cos(#) axis. Bottom right: projection of the "upstream minus 
downstream" distribution. Solid markers are for points inside the fiducial region, 
empty markers for outside. The contours delimit the fiducial region. See text. 
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Figure 9.2: Fit residuals for set B. Top left: color-coded residuals in the cos(0) 
vs momentum plane. Top right: projection on the momentum axis. Bottom left: 
projection on the cos(#) axis. Bottom right: projection of the "upstream minus 
downstream" distribution. Solid markers are for points inside the fiducial region, 
empty markers for outside. The contours delimit the fiducial region. See text. 
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Figure 9.3: Fit residuals for set 1.96 T. Top left: color-coded residuals in the cos(#) 
vs momentum plane. Top right: projection on the momentum axis. Bottom left: 
projection on the cos(#) axis. Bottom right: projection of the "upstream minus 
downstream" distribution. Solid markers are for points inside the fiducial region, 
empty markers for outside. The contours delimit the fiducial region. See text. 
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Figure 9.4: Fit residuals for set 2.04 T. Top left: color-coded residuals in the cos(#) 
vs momentum plane. Top right: projection on the momentum axis. Bottom left: 
projection on the cos(#) axis. Bottom right: projection of the "upstream minus 
downstream" distribution. Solid markers are for points inside the fiducial region, 
empty markers for outside. The contours delimit the fiducial region. See text. 
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Figure 9.5: Fit residuals for cloud muon set. Top left: color-coded residuals in the 
cos(#) vs momentum plane. Top right: projection on the momentum axis. Bottom 
left: projection on the cos(#) axis. Bottom right: projection of the "upstream minus 
downstream" distribution. Solid markers are for points inside the fiducial region, 
empty markers for outside. The contours delimit the fiducial region. See text. 
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Chapter 10 

Conclusion 



The TWIST result for 5, Eq. (|9.1|) . is consistent with the Standard Model prediction 
5 = | . It is also consistent with the previous best measurement [27[ E3 5 = 
0.7486 ± 0.0026 (stat.) ± 0.0028 (syst.) . TWIST result Eq. flST} can be rewritten 
with the errors combined: 

5 = 0.74964 ± 0.00130. (10.1) 

Compared to the combined error of [27], the TWIST result is an improvement of a 
factor of 2.9. Because the measured value is consistent with the Standard Model, it 
places new limits on possible deviations from the theory. 

10.1 Model-independent limits 

on right-handed muon interactions 

Model-independent limits on right-handed couplings of the muon can be obtained 
using Eq. (flOTj) . the TWIST result [98] 

p = 0.75080 ± 0.00105, (10.2) 

and a value of P^S/p. Using [TOOl [TOU HD2J we get 1 

P^S/p = 0.99787 ± 0.00082. (10.3) 

No erratum correcting a mistake in p, — e scattering 100 has been published for [101] . In 
[T02] . page 103, the value of P^S/p = 0.9984 ± 0.0016 ± 0.0016 is quoted. Removing an upward 
correction factor of 1.0007 (page 86) for depolarization in p — e scattering, we obtain P^^S/p — 
0.9977 ±0.0016 ±0.0016. The latter number, combined with P^5/p = 0.99790 ± 0.00046 ± 0.00075 
from [100] . gives the value quoted in the text. Because [100111011 1102] quote their final results not 
as values but only as lower limits on P^^S/p, Eq. (]10.3[1 can only be used to produce limits on, but 
not values of, other parameters. 
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We can transform Eq. ()1.14j) 

<?H{i + 5e-irtf} (10 ' 4) 




(10.5) 



Since 6 » 3/4, 16 5/9 - 1/3 > 0. Also P M < 1, therefore 

<&s\{i-K(T'-k)} < 10 - 6 > 
-5{i-W//>)p("-s)}- <ia7) 

Substituting the values of P^S/p, p, and 5 

Q R = 0.00061 ± 0.00086. (10.8) 

Because mathematically Q R > 0, we convert Eq. (|10.8|) to a 1-sided limit: 

Q R < 0.00184, 90% confidence level. (10.9) 

This is our new limit on the fraction of muons decaying through right-handed inter- 
actions. 

Using Eqs. (|1.10p - (|1.13jl and (|10.9p . we can put new limits on interactions that 
couple right-handed muons to left-handed electrons. These limits are summarized 
in Table [TaTI 

Coupling TWIST limit Previous limit from [87] 
\g s LR \ 0.086 0.125 

\g v LR \ 0.043 0.060 

\gl R \ 0.025 0.036 



Table 10.1: 90% confidence level upper limits on couplings between right-handed 
muons and left-handed electrons. 

10.2 Limits on 

From the same inputs of 5, p, and P„£6/p, as in the previous section, it is possible to 
place new limits on Pu£. Using fllQ.3[) . fllQ.2[) . and (jlO.lj) . we obtain an intermediate 
result in the form P M £ = vl ± 0L, which is not a "value", but can be used to set 
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a lower limit on (see footnote on page I80[) . Upper limits on £ are imposed by 
£<Vp — 1 an d > 0, with the latter being the strongest. From Q R = we get 
£ = vjj ± ov, and can use it to compute a limit. Because P M £ < £, the final result 
should be in the form L < P^£ < £ < [/. This makes L and {7 weaker than the 
corresponding one-sided limits would be. We defined the lower, L, and the upper, 
U, bounds, as L = vl — kaL, U = vu + kajj, imposing the same sigma multipliers 
at both ends. Demanding that the sum of integrals of the normal distributions 
Gauss(t>L, ol) and Gauss (vjj, ay) between L and U is 2 x 0.9, we obtain a "90% 
confidence interval" 

0.9960 < P^ < £ < 1.0040. (10.10) 

(In fact the two integrals are 0.8975 (L), and 0.9025 (U), so that each of the limits 
is very close to 90%.) 

10.3 Limits on left-right symmetric models 

The lower limit on P^£ (jlO.lOp can be used to put limits on mass of the second 
charged gauge boson and its mixing angle with the Standard Model W in left-right 
symmetric theories, see (|1.2ip - (jl.22p . 

Manifest left-right symmetric models assume that gi = gR, the right-handed 
CKM matrix coincides with the known left-handed CKM matrix, and there is no 
CP violation in the mixing: uj = 0. (Notation from section fl.3.21 ) Pseudo-manifest 
models allow CP violation, but require V L = (V R )* and gi = gR. See pp. 377-379 
in [87] for a recent review. 

Some of the existing constraints in the mass — mixing angle plane for the manifest 
case are shown on Fig. 110. li It can be seen that TWIST constraints indeed provide 
an improvement over previous muon decay data, including a dedicated search [100] . 
The new limit on Wr mass is m,2 > 420 GeV/c 2 , compared with the previous limit of 
406 GeV/c 2 [100] (402 GeV/c 2 with the modern value mi = 80.423 GeV/c 2 ). This 
new limit is also significantly stronger than the combined limit from many nuclear 
beta decay experiments summarized in [103] . 

A measurement of the Michel parameter p provides a constraint on the mixing 
angle, which does not depend on the mass. The best limit, established by TWIST, 
is I C| < 0.030 at 90% confidence level [98]. There is also a very tight constraint 
on mixing angle from superallowed nuclear beta decays [105] . that is dependent on 
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Figure 10.1: 90% CL constraints on manifest left-right symmetric models. The 
excluded region is low mass and large |£|. Bold line: TWIST lower limit on P^. 
Dashed line: dedicated search for right-handed currents in muon decay [lOQj . Dash- 
dotted line: one-sided limit P^ > 0.991867 (90% CL) from a direct measurement 
of P M £ |104| . Dotted line: combined nuclear beta decay data |103j . 

nuclear theory and other inputs. When PDG recommended values [87) are used for 
elements of the CKM matrix, it yields a non-zero result ( = 0.00176 ± 0.00074. 

Limits from direct searches at colliders do not constrain but provide better 
constraints on the mass. The strongest combined result is ni2 > 786 GeV/c 2 at 
95% confidence level [HI]. (At 95% CL, 0.99523 < P^ < £ < 1.00472, and m 2 > 
402 GeV/c 2 from TWIST.) Collider results need less restrictive assumptions about 
mass of right-handed neutrino than low-energy tests, but depend on the assumed 
decay channels of the right-handed boson. 

Under the assumption of manifest left-right symmetry a yet stronger limit of 
m2 > 1-6 TeV/c 2 can be extracted from the Kl — K$ mass difference Amjf [107] . 

Manifest and pseudo-manifest left-right symmetric models have severe difficulties 
explaining experimental data (see e.g. j!09j ). therefore a more general case has to be 
considered. The number of parameters become much larger in generalized models, 
since no statement about the right-handed CKM matrix can be made. Many limits 
become significantly weaker for generalized left-right symmetric models, and may 
cease to be useful if fine-tuning is allowed. This is true for constraints from Attik 
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Figure 10.2: A comparison of model dependence of muon decay and collider limits 
on left-right symmetric models. Solid lines: 95% CL limits on manifest left-right 
symmetry from TWIST (curved line) and DO [108 J (vertical). Dashed lines: TWIST 
exclusion for arbitrary CP violation and right-handed CKM matrices (curved), DO 
mass limit for a specific set of parameters of non-manifest model (vertical). 

[1101 1109] • nuclear beta decays (Fig. 12 in [103] ). and collider results [lllj . Muon 
decay data, on the other hand, have very little sensitivity to assumptions about the 
unknown right-handed sector. An illustration of this statement is given on Fig. 110.21 
The two vertical lines represent 95% DO limit under assumption of manifest sym- 
metry (solid), and a significantly weaker one obtained assuming a different specific 
set of model parameters (dashed) [108] . To make a direct comparison to the quoted 
DO results, the TWIST limit from Fig. 110.11 was converted to 95% CL, and shown as 
the solid curved line. The dashed curved line on Fig. 110.21 shows 95% CL excluded 
region from TWIST for the most general case, with an arbitrary fine tuning allowed. 
If a point is outside of the dashed curve, it is excluded for any possible combina- 
tion of the parameters with at least 95% confidence level. (Of course, right-handed 
neutrinos still have to be light.) 

A very strong limit on the mass of the second W boson, of the order of 3 TeV/c 2 , 
can be obtained from big bang nucleosynthesis [112] , They require right-handed 
neutrinos to be lighter than about 1 MeV/c 2 , and depend on assumptions about 
cosmological models that are outside of the scope of particle physics. 
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Among tests of left-right symmetric models done within the context of particle 
physics, muon decay provides better mass limits than nuclear (and neutron) de- 
cay data. These limits are still weaker than limits from direct collider searches. 
Unlike direct collider searches, muon decay constrains the mixing angle as well as 
the mass parameter. But muon decay constraints are conditional on the lightness 
of right-handed neutrinos. On the other hand, collider results are obtained using 
assumptions of manifest or pseudo-manifest symmetry. An important advantage of 
muon decay results is their weak dependence on unknown parameters of more gen- 
eral left-right symmetric models, making them complementary to data from other 
sources. A discussion on complementarity of different observables for generalized 
left-right symmetric models can be found in [103] . Muon decay results are also not 
subject to complications of QCD and nuclear theory. 

10.4 Limit on non-local tensor interactions 

Using (|1.24p and a 90% confidence level lower limit for 5 from (jlO.ip . we obtain 

\g% R \ < 0.024, 90% confidence level. (10.11) 

Since the proposed value is g RR « 0.013 (Section ll.3.3p . this limit does not constrain 
the model significantly. 
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Appendix A 

TWIST coordinate system 
and kinematic variables 

The origin of the TWIST coordinate system is in the center of the muon stopping 
target. The Z axis is along the detector stack, in the direction of the muon beam 
(see Fig. 12. 3p . The Y axis points upwards, and the X axis completes a right-handed 
XYZ coordinate system. The U and V axes are obtained from the X and Y 
axes, correspondingly, by a +45° rotation about the Z axis. TWIST wire chambers 
measure U and V coordinates, not X and Y, see chapter [2j 

The angle 6 of a track is defined by cos(#) = p z /p, where p = \p\ is the momentum 
of the particle, and p z is the projection of the momentum on the Z axes. The 
transverse momentum, pt, is defined as pi = p 2 — p\. 

The "upstream" and "downstream" directions are defined relative to the muon 
beam, that is cos(#) < for an upstream decay, cos((9) > for a downstream decay. 
Upstream part of the detector is the one seen by an incoming muon before it comes 
to rest in the central stopping target. 
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Appendix B 

Optimization of fit range 
for the energy calibration 
procedure 



The energy calibration procedure attempts to compensate for differences between 
data and Monte-Carlo that affect the position of the end point. To accomplish that 
task, the calibration results should not be sensitive to these same differences. This 
can be re-stated as a requirement that the fit should be able to recoup a change in 
f3, a u , ad- So the optimization criteria is not the minimization of a fit bias, but the 
minimization of any dependence of a bias on the shape of the end point region. 

To quantify the ability of the fit to recoup a shape change, a data spectrum 
was distorted by applying the energy calibration transformation (I6.10p with e.g. 
P t = 25 keV/c. The energy calibration procedure was run on both the original 
and the distorted spectrum, and the changes 
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were computed. These changes characterize stabilities of different parameters to 
the given shape change. To obtain more reliable estimates, several values of /3 shlft , 
from —75 keV/c to +75 keV/c in steps of 25 keV/c, were used, and RMSes of A/3, 
Aa u , Aa d , Aa computed. RMSes of Aa sum and Aa^is were also obtained from 
the same data. Similar scans were done for a u and a^, to quantify stabilities of 
fit results under different shape changes. A "variation" of ao was accomplished by 
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smearing the reconstructed momentum with a Gaussian, whose width was defined 
following (|6.7p aso-g hift /|sin(6»)|. For that scan, Ao"o was defined through a quadratic 
difference Aa = y /, (<rg hifted ) 2 - (cr§ hift ) 2 - a 1 ^. Therefore, for a fixed choice of fit 
range, we had 24 numbers, characterizing stabilities of 6 fit results (/?, a u , a^, (Tq, 
CKsum, CKdiff) under 4 different shape changes (the scans of /3 shlft , a^ hlft , a| hlft , <7o hlft ). 
Three versions of choosing the momentum range were tested: 

Pedge(fl) - Ci < p < Pcdgo(^) + C 2 , (B.5) 
Pedge(0) - Ci < p < PcdgcW + «2C r (6 l ), (B.6) 

PcdgcW - si<r(6) <p< Pcdgc(^) + s 2 a(9). (B.7) 

where Cj are constant momentum intervals and are constant multipliers. 

For each of the schemes (|B.5|) - (|B.7P a 2-dimensional scan of the parameters q 
and/or Sj was performed, with c\ = 0...2.5 MeV/c in 0.25 MeV/c steps, C2 = 
... 0.5 MeV/c in 0.05 MeV/c steps, si = . . . 5 in steps of 0.5, and s 2 = ... 5 in 
steps of 0.5. At each scan point, the 24 "stability" numbers were computed, so that 
a scan yielded 24 2-dimensional "maps" of the fit range parameter space. 

These 24 maps were examined by eye, and a "compromise" region, approximately 
minimizing all of the stability parameters, was identified. Then the best regions 
found for flB.5[) — fjB.Tj) were compared to each other. The best results were given by 
the scheme (|B.6|) . with c\ = 0.75 MeV/c and s 2 = 0.5. So these were the settings 
used for the energy calibration during production fitting. 
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Appendix C 

The spectrum expansion 



Integrating the response function K from (|7.ip over x in bin i of the spectrum 
histogram, we get a "binned" response function Ki(x'), that is, the probability to 
get the reconstructed event in bin i. For a dataset with N' true decays, the expected 
number of events reconstructed in bin i is: 

Ni(X) =B i + N' J Ki(x') f(x'; A) dx' (C.l) 

where is the background, f(x'\ A) is the true distribution of decays, and SIq is the 
whole kinematically allowed phase space. Often we have an analytical expression 
for the theoretical distribution representing / only up to a normalization factor: 

f(x'; A) = A(X) Fix'; A), J f(x'; A) dx' = 1, (C.2) 

9.0 

we know F but do not know A. (Of course A can be calculated numerically.) This is 
true for muon decay: (jl.3p gives the differential decay rate F, but not the probability 
distribution f. 

TWIST measures the shape of the spectrum. We can get rid of the absolute 
count by normalizing iVj to the total number of reconstructed events in fiducial 
volume f2: 

m(X) = Ni(X)/N(X), where N(X) = JV 4 (A). (C.3) 

n 

A change in the parameters A modifies the spectrum shape as 

n t (X + AA) - n,(A) = AA 1^ - n,(A) AA ^ 1^ + 0(AA 2 ) (C.4) 

n 
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Using (|C.1|) - (|C.3|) we can express the derivative on the right hand side as: 
AX 1 dNi N' . f v . ., dF(x';X) , , 

Ni A , 1 3i4 A x 1 dA l2 



(C.5) 



+ ¥ AA I(A)aA-AF AA I(A)9A + ^ AA 



Also, 



no 

and 

A" 1 (A) = / F(x';X)dx'. (C.7) 



n 

Analytical expressions for F and dF/dX are known, so there are available many 
ways to calculate integrals (|CJ.6[) — (|(J.Tj> . The integral in (|C.5p contains an unknown 
response function K, and therefore must be evaluated with Monte-Carlo. It is 
convenient to calculate (|C,6p - (jC,7p by Monte-Carlo integration as well. 

Doing the integrals 

A definite integral of a bounded non-negative function g(x') can be evaluated us- 
ing the acceptance-rejection method: choose a y max > max x'en 9( x ') an d sample 
-^thrown points {x' , y} from a uniform distribution on Qq x [0,Y max ]- Call a point 
accepted if y < g(x'). Then 

g(x') dx' » Y max / dx'. (C.8) 



Athrown 

where A acc is the count of accepted points. This recipe is applicable for evaluat- 
ing (|C.7p . An integral for a more general g{x') can be written 

g(x')dx' = J g + (x')dx' - J g~(x')dx' (C.9) 

Q,q S~2q 17q 

where <? + and g~ are non-negative functions: 

JsOO ifg(x')>0, _ f\g(x')\ i{g(x')<0, 

g [x ) = < 3 (a; ) = < (C.IO) 

otherwise, otherwise. 
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Then from (l08l) 



g(x') dx' « ^% ^ y max / dx', 

1 ^thrown 



where 



F max > max (C.12) 
y<g + (x>) for iV+ c , (C.13) 
y<g-{x') for iV- c . (C.14) 

The integral in (|C.6p can be calculated in this way. 

To evaluate the integral in (|C.5P we can use the following procedure: with the 
"generating function" g{x') = dF{x';\)/d\ sample {x',y} from a uniform distri- 
bution as before. For every point accepted according to (|C,13P or (|C.14p use x' to 
define the kinematics of an event in GEANT, simulate and reconstruct the event. If 
the event passed the whole analysis chain and landed in bin i, count it as N^[^] 

d- Then 

/ K iW) ™ « = N >f ~ ^Igi y_ [fSl / * (C, 5 , 

J OA -^thrown 77Y ■/ 



or iVr[f ]. Then 



and 



^thrown \rsx J 

120 "0 



A ,\ 1 ^ = A A ymax[ ^ ] thrown [i 7 ] 

TV y max [F] N[F] 



{ AWowntlf] ^ ^thrown [§f] J 

+ 0(AA 2 ). (C.16) 

N' here, as before, denotes the number of "true" decays, that is, those accepted 
in generation but not necessarily passed through the analysis. The argument in 
the square brackets indicates which generating function was involved. For example, 
Ni[F] = N u and N'[F] = N' from flgHp . The numbers N, N u iVf , N' , 7V /=fc are to 
be understood as statistical expectations. 

Substituting (IC.16P into (IC4[) and introducing 

m^] = Ntm - jvrffl, jv'i^ ] = ^[f] - ansa, (c.i7) 
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we finally get 

m (X + AA) - m(X) = AA Z^Iffi thrown [F] 
{ } K ' Y max [F] N[F] 

\iV thrown [ff] iV'[F] iV thrown [|f] iV[F] iV thrown [|f] 



iV[F] iV'[F] iV thrown [§f] 
where B = ^2 ie ^ Bi. Introducing the efficiency 



+ 0(AA 2 ), (C.18) 



'W-^lCT <c ' 19) 



of the mapping of Qq x [0, Y max ] into Q, the integral of dF/dX 



N'\—] 

v{x) = y max [f ] ld y dF , (c.20) 



and normalized spectra 



OX J Ar _ r dF 

we can re-write (|0. 18|) as 



Vi(\)=Y max m] ,L ^_ , ft = ^, (C.21) 



nj(A + AA) 



nj(A) 



1- J^AA a £:- 1 (^ Q + /5^ a ) 

m 

+ ^ AA a V? + &Z> a ) + 0(AA 2 ). (C.22) 

a=l 

Here the index a = 1, . . . , m numbering the components of A is shown explicitly, 

P = Eign ft. and ^ a = Ei 6 n 

Equation (|C.22p shows that a reconstructed spectrum for parameter values A + 
AA can be represented as a linear combination of reconstructed spectra with different 
values of parameters A. The coefficients in front of the "derivative" terms in the 
linear combination are proportional to the deviations of parameters AA. 
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Appendix D 



Conversion formulas for the 
P^S parametrization 



Here are the formulas for conversion between the Pu£, 5, covariance matrix V, and 
z = Pfi^\p^s,w = Pfi£,$, covariance matrix U, parametrizations. 

1) P^S^P^S 
Az = AP^ 

Aw = (P^ + AP^)(S + A5) - P^ 5 

Uij = Vij, i,j^w 

U lw = (d + AS)Vit + (P^ + AP^)V iS , i + w 

U ww = (6 + A,5) 2 % + 2(P^o + AP M (6 + A5) V^s + (P^o + AP^f V S5 

2) P^5^P^,S 

AP^ = Az 

A5 = (Aw-5 Az)/(P^o + Az) 

Vij = U^, i,j^S 

P^ 6 + Aw TT 1 ., x 

(P^o + Az) 2 P^ + Az 
_ (P^q5 + Aw) 2 TT P^qSq + Aw 1 
V&S ~ {P^ + Azf Uzz (P M & + Az)* Uzw + (P^ + Azf Uww 

The covariance matrix conversion is approximate, see e.g. [95] for details. 
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